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Phase transitions from an active into an absorbing, inactive state are generically described by the 
critical exponents of directed percolation (DP), with upper critical dimension d c — 4. In the frame- 
work of single-species reaction-diffusion systems, this universality class is realized by the combined 
processes A — > A+A, A + A — > A, and A — > 0. We study a hierarchy of such DP processes for particle 
species A,B, . . ., unidirectionally coupled via the reactions A — ► B, . . . (with rates fiAB, ■ ■ •)• When 
the DP critical points at all levels coincide, multicritical behavior emerges, with density exponents 
Pi which are markedly reduced at each hierarchy level i > 2. This scenario can be understood on 
the basis of the mean-field rate equations, which yield f3i = 1/2 I_1 at the multicritical point. Using 
field-theoretic renormalization group techniques in d — 4 — e dimensions, we identify a new crossover 
exponent <j>, and compute cj> — 1 + 0(e 2 ) in the multicritical regime (for small ^ab) of the second 
hierarchy level. In the active phase, we calculate the fluctuation correction to the density exponent 
on the second hierarchy level, = 1/2 — e/8 + 0(e 2 ). Outside the multicritial region, we discuss 
the crossover to ordinary DP behavior, with the density exponent (3i = 1 — e/6 + 0(e 2 ). Monte 
Carlo simulations are then employed to confirm the crossover scenario, and to determine the values 
for the new scaling exponents in dimensions d < 3, including the critical initial slip exponent. Our 
theory is connected to specific classes of growth processes and to certain cellular automata, and the 
above ideas are also applied to unidirectionally coupled pair annihilation processes. We also discuss 
some technical as well as conceptual problems of the loop expansion, and suggest some possible 
interpretations of these difficulties. 

PACS numbers: 64.60. Ak, 05.40. +j, 82.20.-w. 



I. INTRODUCTION 

The notion of universality plays a central role in equi- 
librium as well as in non-equilibrium statistical mechan- 
ics. It was first used by experimental physicists in order 
to describe the observation that certain thermodynamic 
observables measured in different and apparently unre- 
lated equilibrium systems near a continuous phase tran- 
sition may exhibit the same type of singular behavior jl| . 
It was, in fact, then realized that the majority of equi- 
librium critical phenomena belong to very few univer- 
sality classes which are characterized by a certain set of 
critical exponents. In order to explain universality, vari- 
ous theoretical approaches have been constructed, for ex- 
ample scale invariance Q , field-theoretic renormalization 
group techniques , and the theory of conformal invari- 
ance j|] , which predicts a series of universality classes for 
two-dimensional critical systems. Thus in equilibrium 
statistical mechanics, especially in two dimensions, the 
concept of universality seems to be well understood. 

For systems far from equilibrium, however, the situa- 
tion near dynamic continuous phase transitions is less 
clear. Non-equilibrium processes are much harder to 
solve or even to characterize exactly since the probabil- 
ity distribution cannot be obtained from an energy func- 
tional, but has to be derived directly from the equations 
of motion. In addition, systems far from equilibrium are 



in general not conformally invariant since there is no sym- 
metry between spatial and temporal degrees of freedom. 
Nevertheless it appears that universality, although prob- 
ably in a weaker sense, may also play an important role 
in non-equilibrium critical phenomena. As in the case of 
equilibrium physics the picture that emerges is that only 
a few distinct universality classes seem to exist. The 
known examples include phase transitions in driven dif- 
fusive systems J5j, the power-law decay in annihilation- 
coagulation processes the "parity-conserving" dy- 
namic transition for branching and annihilating random 
walks with even offspring number ||, non-equilibrium 
roughening transitions in growth models [ fl0[ , specifically 
in the KPZ equatio n |llfl , and the critical points of di- 
rected percolation ]12(|, as described by Reggeon field 
theory |l3), and of dynamic (isotropic) percolation p4| . 
As a unifying theoretical framework is not yet available, 
we are still far from a systematic classification of non- 
equilibrium critical phenomena. Therefore one impor- 
tant direction of research is in fact to search for further 
unknown universality classes. 

Another direction, which is actually the objective of 
the present paper, would be to investigate the known 
universality classes in more complicated contexts. The 
basic idea is to use several non-equilibrium systems of a 
known universality class as building blocks of a superior 
structure in which the systems are linked to each other in 
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a specific manner. The question posed is whether these 
systems can be combined in such a way that novel critical 
behavior emerges. In other words, is it possible to couple 
several non-equilibrium systems of a given universality 
class in such a manner that the resulting critical behavior 
is characterized by independent critical exponents? 

Such a model with quadratically coupled directed per- 
colation (DP) processes was recently investigated by 
Janssen, who found that despite the apparent complex- 
ity of this coupled multi-species system, the universality 
class of the active / absorbing transition was that of DP 
itself |l5| . In the present article we show that novel crit- 
ical behavior may, however, occur when several copies of 
the same non-equilibrium process are linearly coupled in 
one direction without feedback. More precisely, we con- 
sider a linear hierarchy of unidirectionally coupled copies 
A,B,C, ... of the same non-equilibrium system: 

B -^C -> ... (1.1) 

The systems are coupled in such a way that the dynam- 
ical processes at a certain level in the hierarchy depend 
on the state at the preceding level but not vice versa. For 
example, subsystem A, the lowest level in the hierarchy, 
is not influenced by the dynamics of B and C and thus 
it evolves independently as if the other hierarchy levels 
did not exist. Subsystem B in turn is affected by A but 
not by C, and hence this is the first level in the hierarchy 
where novel critical behavior might occur. The hierar- 
chy can be continued to infinitely many levels. However, 
because of the unidirectional structure one can always 
truncate the hierarchy at some level without affecting 
the temporal evolution at lower levels. For example, we 
may consider a two-level hierarchy A — > B or a three- 
level hierarchy A — ► B — * C; in both cases the dynamics 
of the subsystems A and B would be exactly the same. 

The most interesting behavior of the composite sys- 
tem is expected when the subsystems A, B, C, . . . them- 
selves are close to criticality. This usually happens when 
the isolated subsystems would undergo a continuous non- 
equilibrium phase transition. Let us assume that the 
stochastic process under consideration is controlled by a 
single parameter p with a phase transition taking place 
at p — p c . A unidirectionally coupled hierarchy of 
such processes is thus controlled by a sequence of in- 
dependent control parameters p^ A \p^ B \p^ c \ . . ., which 
means that the composite system is described by a high- 
dimensional phase diagram. When all levels are critical 
(i.e. p( A ^ — = = . . . = p c ), a complex interplay 
of long-range correlations is expected. We will refer to 
this special point in the phase diagram as a multicriti- 
cal point. In the vicinity of this point the properties of 
the entire system depend crucially on the direction from 
which it is approached, resulting in interesting multicrit- 
ical behavior. In particular we will consider the special 
case p( A ^ = p( B ^ = p( c ' = ... = p, where the entire 
hierarchy is controlled by a single parameter. 

The outlined concept of unidirectionally coupled non- 
equilibrium processes is quite general and may be applied 
to various dynamical systems. But, as we shall also see, 



this mechanism does not necessarily lead to new univer- 
sality classes in all cases, and we have already mentioned 
the quadratically coupled DP processes jl5| as one coun- 
terexample. In the present work we will focus on non- 
equilibrium processes which display a continuous phase 
transition from a fluctuating into an absorbing state, i.e. 
a configuration which once reached, cannot be escaped 
from. 

The canonical example for a transition into an absorb- 
ing state is the critical point of directed percolation (DP). 
In DP, sites of a lattice are either occupied by a particle 
(active) or empty (inactive). The dynamic processes are 
that a particle can self-destruct or produce an offspring 
at a neighboring empty site. If the rate for offspring pro- 
duction p is very low, the system always reaches a state 
without particles which is the absorbing state of the sys- 
tem. On the other hand, when p exceeds a certain crit- 
ical value pc, another steady state with a finite particle 
density exists on the infinite lattice. In between a contin- 
uous phase transition takes place which is characterized 
by long-range power-law correlations. Another example 
considered in the present work is the annihilation process 
A + A — > in which the particle density decays as 
f-d/2 m dimensions d < 2. Here the absorbing state (the 
empty lattice) is approached without the tuning of any 
parameter, i.e. the process is "critical" by itself. 

To construct a unidirectionally coupled hierarchy of 
such processes we implement additional dynamical rules 
which allow each particle at a given level in the hierar- 
chy to induce the creation of a new particle at the same 
lattice site of the next level. This ensures that the com- 
posite system still has an absorbing state, namely the 
empty state without particles. More precisely, a whole 
hierarchy of absorbing subspaces is generated. For ex- 
ample, if subsystem A enters the inactive state, it will 
never become active again and therefore the dynamical 
processes are restricted to an absorbing subspace where 
only B,C, . . . may fluctuate. This subspace in turn con- 
tains another absorbing subspace in which B is inactive, 
etc. As we will demonstrate, such a hierarchy of sys- 
tems with absorbing states coupled by induced particle 
creation is characterized by a subset of novel critical ex- 
ponents. It should be emphasized that the emergence 
of novel critical behavior is related to the fact that the 
processes are coupled in only one direction. Even a very 
small feedback (e.g. B — > A, C — > B) or cyclic closure 
(e.g. A — > B — > C — > A) would destroy this new feature. 

In this article we present a detailed analysis of unidi- 
rectionally coupled DP [M . For the case of equal con- 
trol parameters it is observed that the asymptotic par- 
ticle densities near the multicritical point are character- 
ized by different critical exponents Pa, Pb, Pc, ■ ■ ■■ Since 
level A evolves independently, Pa is just the usual den- 
sity exponent of DP. At higher levels numerical estimates 
show that the density exponents are considerably reduced 
compared to their DP values. In Sec. II we discuss the 
mean-field theory of coupled DP which already explains 
why the density exponents at higher levels are reduced. 
It also allows us to study crossover phenomena close to 
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the multicritical point. Mean-field theory is expected to 
hold above the DP critical spatial dimension d c = 4. 
For d < d c the numerically observed values for the den- 
sity exponents are much smaller, which means that the 
mean-field results are strongly modified by fluctuation 
effects. In order to understand these reduced values, we 
recently derived the critical exponents to one-loop or- 
der |L6| by means of a field-theoretical renormalization 
group analysis, based on the "Hamiltonian" representa- 
tion of the classical master equation p jl7| , |i^| . In Sec. Ill 
we present these calculations in detail, including a dis- 
cussion of the diagonalized and multicritical theories, re- 
spectively, the active phase, logarithmic corrections at 
d c = 4, crossover studies, and the critical behavior at 
higher levels in the hierarchy. The field-theoretical re- 
sults are supported by extensive numerical simulations in 
Sec. IV, while in Sec. V various applications of coupled 
DP are demonstrated. The examples of coupled annihi- 
lation and other closely related topics are the subjects of 
Sec. VI. Finally, a critical discussion of some technical as 
well as conceptual problems arising in the field-theoretic 
approach are the subject of our conclusions in Sec. VII. 



II. COUPLED DP PROCESSES: MEAN-FIELD 
APPROXIMATION 

In the bulk of this paper, we shall study unidirection- 
ally coupled directed-percolation processes. It is conve- 
nient to represent these processes in the framework of 
reaction-diffusion systems. The starting point of the hi- 
erarchy of coupled systems is therefore the following re- 
action scheme, which can be viewed as the prototype for 
the active / absorbing state transitions in the directed- 
percolation universality class [p~3| : 



A — > A + A with rate a A 
A — ► with rate /i A 
A + A—> A with rate X A 



(2.1) 
(2.2) 
(2.3) 



We can immediately write down the corresponding rate 
equations for the particle density n A (t); this description 
neglects fluctuation and correlation effects, and therefore 
corresponds to a mean-field approximation. The average 
part icle number is increased via the b ran ching reaction 
(2.1), and reduced via both the decay ( |2.2D and coagula- 
tion (2.3). However, while the first two processes occur 
spontaneously, with rates a a and /j. a , respectively, and 
therefore the change in particle number is proportional 
to the particle density itself, the coagulation reaction re- 
quires that two particles A meet on the same lattice site 
(in a discrete representation) , an d hence the total par- 
ticle loss due to the process ( |2.3|) is proportional to the 
density squared. This yields the balance equation 



dn A (t) 
dt 



(a a - HA)nA{t) - X A n A (tf 



(2.4) 



As we are considering local reactions only, we may gen- 
eralize this mean-field equation slightly by considering a 



coarse-graine d loc al particle density n A (x, t), and supple- 
menting Eq. (2.4) with a diffusion term, 



dnA ^ l) =D(V 2 - r A ) n A (x, t) - X A n A {x, tf . (2.5) 

Here we have introduced the diffusion constant D, and 
defined r A = (ft a — &a)/D. 

Obviously, the mean-field dynamic phase transition oc- 
curs at the point = 0, where the balance of gain and 
loss due to the processes linear in ua ch ange s sign. For 
rA > 0, the only stationary state of Eq. (2.5) is tia = 0, 
and for t — > oo the particle density will simply decay to 



zero according to n A (t) — » e 



-Dr A t 



because once the par- 



ticle density has become sufficiently small, the coagula- 
tion contribution can be neglected. Furthermore, n A = 
represents an absorbing phase, because once there are n o 
particles left in the system, none of the processes (2.1)- 
( |2.3D can happen any longer — hence all fluctuations 
cease, and the system cannot escape from this state. For 
ta < 0, on the other hand, there is another stationary 
state with non-zero particle density 



n A = D\r A \/X A 



(2.6) 



which can be viewed as the order parameter of the ac- 
tive phase. In the active state, the asymptotic density 
is approached exponentially again, with the characteris- 
tic rate D|r,i|. Precisely at the transi tion , only the term 
proportional to n A survives in Eq. ( |2.5| ), which there- 
fore becomes identical to the mean-field rate equation 
for diffusion-lim ited coagulation or annihilation Q. The 
solution to Eq. (|2.4|) then becomes 



n A (t) ~ 1/t 



(2.7) 



i.e., the density decays according to a power law at the 
critical point. 

Upon identifying rA = p c —p> where p c denotes the per- 
colation threshold, we can translate the above mean-field 
results to the notation of directed percolation. Eq. (2.5) 
implies that the characteristic length scale £j_ = |rvi| _i/:2 
diverges in the vicinity of the transition, 



U ~ \t a \- 



u x = 1/2 



(2.8) 



At the critical point, the exponential decay rates vanish, 
and the characteristic frequency becomes diffusive, u c ~ 
Dq 2 ; hence 



(2.9) 



Equivalently, upon approaching the critical point, the 
characteristic time scale diverges according to 

£,l ~ \r A \- v t , v\\ =zu ± = l. (2.10) 

Also, at p — Pc, 

n A (t) - t' a , a = 1 . (2.11) 

Finally, in the active phase near the transition, the order 
parameter grows as 
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n A 







(2.12) 



Notice that the exponents a and j3 are related, provided 
the following scaling relation holds, 



n A (r A ,x,t) = \r A fh A (x/£±,t/^) . 

For then at p — p c , h A (0, y) ~ y~Pl v \\ in the limit y 
in order for the jr^j-dependence to cancel, and thus 



(2.13) 



oo, 



0: 



vna = zv±a 



(2.14) 



The above relations therefore define three independent 
critical exponents. Alternatively, the third independent 
exponent may be swapped for r/±, which characterizes 
how the equal-time pair correlation function decays at 
the critical point r A = 0, G(|x|) oc l/|x| rf + z - 2+ '^. 

While the scaling relation ( 2.14 ) remains valid also be- 
low the upper critical dimension, which for DP turns out 
to be d c = 4 (see Sec. HI), the exponent values will be- 
come modified for d < d c as a consequence of strong fluc- 
tuation effects. Directly at the critical dimension, one 
expects logarithmic corrections to the above mean-field 
results. 

The idea is now to combine several (say a number k) 
of such DP processes, i.e., to consider the additional re- 
actions 



B — > B + B with rate erg , 
B — > with rate \xb , 
B + B^B with rate As, 



(2.15) 
(2.16) 
(2.17) 



for particle species B,C, . . ., which are coupled unidirec- 
tionally via the transformation reactions 



A — ► B with rate p, A B , 
A — > C with rate fi A c , 
B — ► C with rate [Ibc > 



(2.18) 
(2.19) 
(2.20) 



etc., but without feedback, i.e., we do not allow processes 
of the type B — > A. This prescription therefore defines 
a hierarchical structure of coupled DP processes. For 
simplicity, we choose identical diffusion constants D on 
each hierarchy level. 

Without the coupling reactions ( [2.1<j ), . . ., for each 
species of particles there is a continuous DP transition 
at = 0, i = A, B, . . .. But the situation changes in 
an interesting way when the transformation processes are 
switched on (except for species A, which is not influenced 
by what happens on the higher hierarchy levels). Let us 
first consider the simplest case of two particle species 
(k = 2) — the generalization to further hierarchy levels 
will t hen be straightforward. The mean-field rate equa- 
tion (2.5) for species A remains unchanged, albeit with a 
modified parameter 



r A 



(pA + f-AB - 0- A )/D 



(2.21) 



Notice that for the first hierarchy level, the sole effect of 
the transformation reactions is an increase of the total 
decay rate ^ = Ma + X); JUi- 



The rate equation for species B, however, contains a 
new gain ter m de scribing the feeding-in of particles via 
the reaction ( |2.1£| ), proportional to the density of A par- 
ticles present. Thus one obtains 



dn B {x, t) 

Ft 



where 



= D(V 2 -r B )n B (x,t) 
-\ B n B {x,t) 2 + ^ AB n A (x 7 t) , (2.22) 

r B = (/is - tr B )/D . (2.23) 



Notice that within the mean-field approximation, n A 
plainly acts as an external source term. Once fluctua- 
tions are important, however, i.e., for d < d c = 4, such a 
simple picture breaks down, especially at the multicriti- 
cal point to be discussed below, where the averages and 
correlations of both n A (x, t) and ns{x, t) are governed by 
power laws. 

We may n ow again search for a stationary solution ns 
of Eq. (2.22), as a function of the mean-field density n A . 
The general solution of the ensuing quadratic equation is 



n B 



Dri 
2Xe 



+ 7 — n A 



1/2 



Dr L 
2\e 



(2.24) 



Thus, for r A > 0, where n A = 0, one finds ub — 
D{\rB\ — ^s)/2As, which is zero for rs > 0, and becomes 
equal to tib = D\tb\/Xb for r B < 0. When species A is 
in the inactive phase, the A and B hierarchy levels are ef- 
fectively decoupled, and we therefore expect an ordinary 
DP active / absorbing transition for species B at rs = 0, 
as in the case \i A b = 0. 



For r A < 0, on the other hand, we have to insert (2.6) 
into Eq. ( 2.24 ). One may now distinguish two situations: 
(i) For (Dr B /2X B ) 2 > D\r A \(i AB /\ A \B, we can approx- 
imate 



n B 



D\r B \ 
2Ar 



D\r A \fiAB ( 2\ B Y 
2X A X B \Dr B ) 



Dri 
2X~e 



(2.25) 



and consequently for rs < 0, we find that ns > 0, and 
the B species is in its active state — the DP transition 
at the half line (r A < 0, tb = 0) present in the uncoupled 
system has disappeared (see Fig. |l|). Instea d, for rs > 
the terms proportional to rs cancel in Eq. (2.25), and 



n B 



\r A \/J, A B/X A r B 



(2.26) 



i.e., the density of species B vanishes as the critical point 
r A = of species A is approached, and with the mean- 
field DP exponent (3=1. Effectively, the DP critical half 
line (r A < Q,rs — 0) for species B has been rotated to 
(r A = 0, rs > 0) in the coupled system. The location of 
the DP critical lines for both species A and B is shown in 
the phase diagram of Fig. [|, where the dotted parabola 
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A active 
B active 



A/B 



A inactive 
B inactive 



A 



A active 
B active 



B 



A inactive 
B active 



FIG. 1. Mean-field phase diagram for the two-level coupled 
DP process. The arrows mark DP active/absorbing transi- 
tions for the A and B particle species, respectively. The dot- 
ted parabola denotes the boundary of the multicritical regime, 
which includes \ta\ = Irs I = \r\ — > 0. 



represents the boundary curve separating the two differ- 
ent regimes for ta < 0. Notice also that in this case 
the " mas s" terms (the contributions linear in ha, n_e) in 
Eq. ( |2.22| ) vanish, and one expects both spatial and tem- 
poral power-law correlations, obviously characterized by 
v± = 1/2 and z = 2. We shall later see that indeed the 
new critical half line {ta = 0, rs > 0) for species B is 
in the DP university class. Summarizing this regime, we 
may say that the B particles are "slaved" by the behavior 
of the A species. 

(ii) The other regime, inside of the dotted parabola 
in Fig. [I], is defined by the condition (Drs/2As) 2 <C 
D\r a\hab I '^a^b ) which includes the special case when 
the critical points of both hierarchy levels are approached 
uniformly, \ta\ — \tb\ = H — > 0. Now we can neglect 
the terms ~ rs in Eq. ( 2.24 ), which yields 



n B 



D\r A \VAB 



1/2 



( Vab 



■ riA 



1/2 



(2.27) 



This implies that the density exponents on each hierarchy 
level i are different in this regime, 



(2.28) 



where 0i = (3a = 1, and 2 — 0b — 1/2 in the mean- 
field approximation. It is to be expected, however, that 
the other independent scaling exponents and z remain 
unaltered; but of course the density decay exponent at 
the critical point Oj wi ll depend on the hierarchy level i, 
according to Eq. (2.14), 



Ui(t) — t 



0i/zvj_ 



(2.29) 



Considering the two-species phase diagram (Fig. [l]) 
again, we see that three critical half lines converge at 
the special point = = 0. The new critical be- 
havior in this regime can therefore be interpreted as the 
effect of a multicritical point 16 1. When this special 



point in the phase diagram is approached along a line 
crossing the dotted parabola in Fig. |l|, one expects a 
crossover from ordinary DP to the new multicritical be- 
havior described by the density exponents 0i and a, . For 
\i"a\ = \tb\ = \r\ — ► 0, the crossover features are all en- 
coded in the generalized scaling function 

n B (r,VAB,x,t) = \rf 1 fiA{\r\~ <l 'iJ,AB/D,x/£x,t/£\\) , 

(2.30) 

where £_l ~ \r\ ~ v± and £ii ~ |r| - "H as in DP. This defines 
the crossover exponent (f>, which constitutes a new scaling 
exponent assoc iated with the coupling [iab- Comparing 
with Eq. (12.24 ); we identify 



1 



(2.31) 



within the mean-field approximation. Furthermore, we 
can use the above mean-field results, which, at the mul- 
ticritical point, imply that n(j/,0,0) ~ y 1 ^ 2 for y — > oo. 
Consequently we have 



02=01- 0/2 



(2.32) 



which relates the new density exponent 02 to t he in dc- 
pendent crossover exponent (j). Of course, Eq. d|3| ) is 
satisfied by the mean-field values. In the fluctuation- 
dominated regime d < d c — 4, however, there will in 
general be 0(e — 4 — d) corrections to both the critical 
exponents and the scaling functions, which will in turn 
lead to a modification of the scaling relation ( 2.32| ). 

One may now readily generalize to higher hierarchy 
levels. For example, for k = 3, one finds the same rate 
equations (2.5) and ( p. 22] ) for species A and B, respec- 
tively, but where now r A = (ft a + Hab + Vac ~ &a)/D, 
and rs — O^b + Vbc — o~b)/D. The mean-field equation 
for nc(x, t) reads 



dnc(x, t) 



dt 

-\cnc{x,t) 2 



= L>(V 2 -r c )n c (x,t) 

I^BcriB(x,t) + HAcnA(x,t) , (2.33) 



which has the general stationary solution 



nc 



(Drc\ 



-,1/2 



VBC 

Ac 



■ nB 



VAC 

Ac 



UA 



Pre 
2Ac" ' 

(2.34) 



A detailed analysis then reveals that, in analogy with 
the two-level hierarchy, there are regions in phase space 
where the C species evolves independently of the lower 
hierarchy levels. On the other hand other regimes exist 
where nc is slaved by either ti^ or ng. In addition, as 
before, a new DP transition may arise for the C particles 
under appropriate conditions, when rA — > or rs — ► 0. 
Furthermore, the previous k = 2 multicritical regime oc- 
curs when either ta > 0, and |t\b| = \ fc \ ~ > simultane- 
ously, or rs > 0, and \ta\ = \rc \ — * 0. 
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As all these features are already contained in our 
above investigation of the two-level coupled DP pro- 
cess, we restrict ourselves to the new behavior emerg- 
ing for k — 3, when all three critical points coincide, 
i.e., \r A \ = \r B \ = \rc\ = M -> 0. More generally, 
in mean-field theory this multicritical regime (with alto- 
gether seven critical quarter planes merging at r = 0) is 
characterized by the conditions r A < 0, (Dtb /2Xb) 2 <C 
D\r A \^ AB /X A X B (as for k = 2), and {Dr c /2\ c ) 2 < 
tiBc(D\r A \fi AB /X A \ B ) 1/2 /Xc + D\r A \fi AC /X A X c . At 
this special po int in par amet er space, n A and ng van- 
ish as in Eqs. (|2.6|) and (2.27), respectively, while 



X a XbXq J 



1/4 



( ^ B C 
V Ac 



n B 



1/2 



(2.35) 



Therefore Ps — Pc — 1/4 in the mean-field approxima- 
tion. Notice also that the indirect A — > C transformation 



rate fi A c does not enter this expression (2.35), which 
only depends on the coupling rates /i A b and /ibc be- 
tween adjacent hierarchy levels. As the latter are not at 
all influenced by the presence of lower levels, this would 
suggest that there exists only one crossover exponent <j> 
describing all the multicritical points generated by the 
unidirectional coupling of the DP processes (with cf> — 1 
in mean-field theory). 

For the density exponents, we conclude that within 
the mean- field approximation on hierarchy level i, using 
2/11 = zv | =1, 



Oi = Pi = 1/2 4 



(2.36) 



For that reason we will employ field-theoretic renormal- 
ization group methods, which allow both for a proper 
derivation of scaling, as well as a systematic e-expansion 
calculation of critical exponents, below the upper critical 
dimension d c = 4. 

Our starting point for a systematic treatment of the 
coupl ed-DP reaction scheme given by ( 2. 1 )-( 2.3) , ( 2.15p - 
( |2.18 ) is an appropriate master equation. On a micro- 
scopic level this comprises an exact description of the 
dynamics. From this equation it is then a straightfor- 
ward process to derive an effective field theory: First 
the master equation is mapped onto a second-quantized 
bosonic operator representation, which is in turn mapped 
onto a bosonic field theory. This procedure is now stan- 
dard, and we refer to Ref. Q for further details. In our 
case, for the two-species coupled DP system we end up 
with the following action 

S = Jd d x Jdt {a [d t + D(r A - V 2 )] a - a A a 2 a + 

+X A (da 2 + a 2 a 2 ) — /i A Bba + (3-1) 
+b [8 t + D(r B ~ V 2 )] b - a B b 2 b + X B (bb 2 + b 2 b 2 )} , 

where we have omitted terms related to the initial state. 
Aside from the taking of the continuum limit, the deriva- 
tion of this action is exact, and in particular no assump- 
tions regarding the precise form of the noise are required. 
Note that if we neglect the terms in the action (3.1) 
quadratic in the response fields a, b, then we recover 
a de script ion identical to the mean-field equations (|2.5|) 



This result should describe the coupled DP multicritical 
point quantitatively correctly for spatial dimensions d > 
d c = 4. Furthermore, on the mean-field level the scaling 
relation ( 2.32| ) generalizes to 



Pi = Pi 



and ( 2.22j ), provided we associate the fields a(x, t), 6(x, i) 
with the coarse-grained local densities of the A, B parti- 
cles. In general, however, below the upper critical dimen- 
sion d c = 4, the terms quadratic in a, 6 (corresponding 
to noise in a Langevin description) cannot be neglected. 

It is now convenient to rescale the fields accord- 
ing to a = (Aa/cta) 1 / 2 ^, a = (a A / X A ) 1/2 if) , b 



I2 l =p x - 0(1 - 1/2 1 ) , (2.37) (\ B /a B ) 1/2 <fo, b = (a B /X B ) 1/2 ^o, and also define new 



compare Eq. ( p. 35] ). Here the observation that only the 
direct transformation rates between neighboring hierar- 
chy levels affect the leading contributio n has entere d cru - 
cially. Agai n, th e mean-field results ( 2.36| ) and (2.31) 
satisfy Eq. ( ^.37| ) tri vially . However, as noted above, 
the scaling relation (2.37) will be modified below the 
upper critical dimension d c = 4, as a consequence of 
0(e = 4 — d) corrections to the scaling function for the 
equations of state, rii (r) . 



III. RENORMALIZATION GROUP 
CALCULATIONS 

A. Preliminaries 

We now turn to a detailed presentation of our field- 
theoretic calculations. As we have pointed out in the 
previous sections, the effects of fluctuations invalidate a 
simple mean-field approach for dimensions d < d c = 4. 



couplings uo = 2(a A X A ) 



1/2 



= 2(a B X B ) 1/2 , as well 



as hq = [i A b{o~ aXb I o~ bXjC) 1 ! 2 (henceforth, the subscript 
"0" denotes unrenormalized quantities). If we introduce 
a length scale k _1 and correspondingly measure times 
in units of k~ 2 (i.e., [Do] — K °)j we find that the new 
fields have scaling dimension K d / 2 , while [r A ] — [tb] — 



[Mo] 



which are thus relevant perturbations in the 



RG sense. On the other hand, [uq] = [u' ] = K 2 ~ d / 2 , and 
the corresponding DP non-linearities become marginal in 
d c — 4 dimensions, as expected |l3|]. It is important to 
note, however, that [X A ] = [Xb] = K 2 ~ d , and hence these 
couplings are irrelevant as compared to uq and u' , and 
may be omitted in the effective action. Finally, we set 
u = uo, such that the theory remains renormalizable 
with equal diffusion constants ]l9| , and consequently ar- 
rive at 



<S'cff = 



d d x I dt {^o [dt + D (r A - V 2 )] iPo 



(V'oV'o - V'oV' 2 ) - AWoV'o + 



"0 



(3.2) 



G 



-ipo [d t + Dq(tb - V 2 )] ip - — (vivo - ¥Wo) } 



We remark that this action is equivalent to the following 
set of coupled Langevin equations 



5 t Vo = A)(V 2 - r A )^ - ^ 2 + ( 



d t tpo = D (S7 2 



rB)fo - y</> + /Wo + V 



(3.3) 
(3.4) 



which represent the obvious a nd e xpecte d gen eralizations 
of the mean- field equations (2.5) and ( 2.22| ), with the 
multiplicative Langevin noise terms 

t)((x',t')) = uoMx, t) S d (x - x')6(t - t') , (3.5) 
(rj(x, t)r)(x',t')) = u ip a (x, t) 8 d {x - x')8{t - if) . (3.6) 

Furthermore, for most of the analysis we will put 
ta = Tb = tq. There are several ways in which the 
effective action fl3.2| ) can be studied. We will begin by 
performing our analysis in the inactive phase, and post- 
pone other methods (including diagonalization and active 
phase computations) until later on. Inspection of the 
above action ( |3.2[ ) reveals that, as expected, the terms 
involving only the ip and i\> fields are exactly the same 
as in the well-known field theory for directed percolation 
(Reggeon field theory) and their renormalization is 
entirely unaffected by the presence of the if and if fields. 
Hence we will begin by briefly reviewing the analysis of 
Reggeon field theory (RFT) in the inactive phase. 



B. DP field theory: Inactive phase calculation 

The renormalization of the RFT action is very well 
known (see Ref. [^0|). The renormalized parameters are 
defined as follows: 



,1/2 

D = Z D D 



ryl/2 J ry -2 

ZJ Ipo T = Z t TqK 



Z u u A d ' k 



1/2 -e/2 



(3.7) 



with e = 4 - d, A d = T(3 - d/2)/2 d - 1 7r d / 2 , and r = 
— foe where tq c is the fluctuation-induced shift of the 
critical point. From the diagrams for the two- and three- 
point vertex functions (see Fig. ^), we can determine the 
one-loop renormalized Z factors. Using dimensional reg- 
ularization and a minimal subtraction scheme, the results 
are 



Zj, = 1 — 



A d K 



Z n = l 



Z u = 1 



8^ 2 e ' 
-»o A d K- e 
IQDl e 
3m 2 , A d K~ e 

16L> 2 e 
hul A d n~ e 

16L> 2 e 



(3.8) 
(3.9) 
(3.10) 
(3.11) 



with roc given by the recursive equation 



r = 



+ 



r 

Vj/ \j/ \|/ 



+ 




+ 




FIG. 2. Diagrams for the two- and three-point vertex func- 
tions to one-loop order for the pure DP field theory. 



roc 



4£> 2 J p r c + P 2 



(3.12) 



where we have used the abbreviation f . . . = 

J P 

J . . . d d p/(2n) d . Defining the flow functions Cu — 
kc* k ln(w/wo) etc., and with an effective coupling v = 
u 2 /16D 2 , the RG (3 function has the form 



V = nd K v = 2u(C„ - Cd) = v(-e + 12v) 



(3.13) 



giving a stable, non-trivial fixed point v* = e/12 + 0(e 2 ). 

The appropriate renormalization group equation for 
the renormalized field (tpR)-, where the angular brackets 
denote averaging with respect to the RFT action, is 



d d d d 

OK 



1 



: C^- 



cV ^ dD dv 2~ 

x(ip R (K,T,D,v,x,t)) = . (3.14) 

Defining the dimensionless field ip as 

{i/)r(k, t, D, v, x, t)) = K d/2 ip(T, v, kx, k 2 Dt) , (3.15) 



the solution of ( 3.14 ), obtained by means of the method 
of characteristics n — > k£, when the coupling v has run 
to its fixed point value v* is 

x^{t£ c - ,v\kx£ 1 k 2 DU 2+c ") . (3.16) 

By inserting the matching conditio n £ — |r| -1 /^^, we 
can now derive the scaling relation (2.13) quoted in the 
previous section. At the fixed point, we obtain 



(if> R (K,T,D,v,x,t)) =\rfi>(v*,- 



KX 



K Dt 



-Vs. ' 



(3.17) 



where the exponents can be identified as combinations of 
the £ functions evaluated at the non-trivial fixed point: 
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12 



+ 0(e 2 ) 



, = JL = 2 + Ci = 2-- + 0(e 2 ) 

/?=^f = ^(d + *-2 + „x) 
= l- e - + 0(e 2 ). 



(3.18) 
(3.19) 

(3.20) 
(3.21) 

(3.22) 



Directly at the upper critical dimension d c — 4 (e = 0) , 
the power laws with these critical exponents are re- 
placed with logarithmic corrections to the mean-field re- 
sults rjx = 0, v±_ = 1/2, — 1, z = 2, and (3 = 1 
(see also Ref. The flow equation for v(£) be- 

comes £dv(i)/d£ = f3 v {£) = 12v(£) 2 , which is solved by 
v(£) — v[l — 12?jln^] _1 , where v — v(£ = 1). Similarly 
£dD(£)/d£ = ( D (£)D(£) = -vD(£)[l - Uvkit]- 1 has 
the solution D(£) = D[l - 12vbx£}^ 12 , and £dr(£)/dl = 
C t {£)t(£) = (-2 + 3v[l - \2v\n£]- l )T{£) is solved by 
t{£) = t£~ 2 [1 - \2v\n£Y x l A . We now integrate the flow 
equations until \t{£)\ = 1, or £ ~ |-r| 1 / 2 (— In l-rD^ 1 / 8 . 
This yields immediately the divergence of both the cor- 
relation length £j_ and the characteristic time scale £y 
upon approaching the phase transition, 



r 2 



rl-V^-inH) 1 / 8 , (3.23) 
/)m- 1 «|r|- 1 (-ln|r|) 1 / 6 . (3.24) 



In order to obtain the corresponding logarithmic correc- 
tion for the density exponent f3, we employ the sol utio n 
(3.16) of the RG equation and the mean-field result (2.6), 
and find 



(3.25) 



where C(£) = exp( f* C,^{£' )d£'/l' ), or equivalently, 
£dC{£)/d£ = £^{t)C(l) = 2vC{£)[l-12v\n£}~ 1 , with the 
solution C{£) = [1 — 12vln_£]_y 6 . Combining everything, 
and setting d = 4 in Eq. ( 3.25| ) finally yields 



|r|(-ln| 



a/3 



(3.26) 



Notice that we had to take care and keep track of the 
dangerous irrelevant variable v here. 



C. Coupled DP field theory: Inactive phase 

We now return to the renormalization of the coupled 
DP field theory. Right from the outset we must take into 
account one key feature of the full theory — namely that, 
on physical grounds, one expects the generation of addi- 
tional mixed cubic vertices. Physically, these novel ver- 
tices correspond to the additionally generated processes 




V 2 
->-< 



-u/2 



-V 2 




FIG. 3. The vertices of the full action S n 



A + B, A^ B + B, A + A^ B, and A + B ^ A, 
with rates gab, v'ab-i ^ab, and X'ab-i sa y- These ver- 
tices must be introduced from the very begin ning , and 
hence we have to replace the above action (3^) with 
•Smc = S c tt + AS, where 



AS 



d d x / dt 



-soipoipo^o ~ y VoVo 
y <Po4>o + s' <f Lp ipo 



(3.27) 



Note that in the shifted theory used above, the new re- 
action processes also modify the bare parameters ta, 
fio, and mq, and furthermore lead to the identification 



•so 



(JAB > 0, 



o~'ab > 0, s ^ — Xab < 0, and 



X'ab — 0- I n the effective Langevin-type descrip- 



tion, Eqs. (3.3)-(3.6) are then replaced by 



dtipo = A)(V 2 



7 2 s , U ,2 , 

d t <p = A)(V 2 - r B )ip - -£-cpl 



so 



where the noise terms satisfy 

C(x,t) = a£i(x,t), 

rj(x, t) = b£ 2 (x,t) + (x, t) 



(3.28) 



(3.29) 



(3.30) 
(3.31) 



Here £i and £2 are uncorrelated white noise variables of 
variance one, and the coefficients a, b and c satisfy: 



a- = MoV'o, 
b 2 = u Q (p - 

c 2 = ^Vc 
u 



^ V'o, 
u 



(3.32) 
(3.33) 

(3.34) 



The complete vertices of the full action S mc are depicted 
in Fig. [|. The propagators for the ip and ip fields are de- 
noted by solid and dashed lines, respectively. The next 
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+ 




+ 



+ — 




+ 




V-->-- + 



y-->- 




FIG. 4. "Mixed" two-point vertex function 
one-loop order for the coupled DP field theory. 



to 



task is to compute the renormalizcd couplings and RG 
fixed points of the full theory. To begi n w ith, we no- 
tice that the DP couplings in the action fl3.2| ) are renor- 
nialized in the same way as in RFT, i.e., the factors 
Zu,, Z T , Zp, and Z u are identical with those in 



Eqs. (3.£ )— ( 3.11 ). Hence we can conclude that the stable 
fixed points for the dimcnsionless rcnormalized coupling 
u = Z u uqA}J 2 n~ e / 2 reads as in DP: 



= [(u/4D)*} 2 = e/l2 + 0(e 2 ) 



(3.35) 



For this reason, the critical exponents rj±, v± and z re- 
main those of the DP universality class for the second 
hierarchy level (i.e., for the B species), and in fact for 
higher hierarchy levels as well. 

However, the same will not be true for the exponents 
(for i > f). For example, the exponent /3% is affected 
by the renormalization of fj,o, an d hence only the expo- 
nent (3\ will remain the same as in DP. In order to com- 
pute the renormalization of we must consider the dia- 
grams renormalizing the "mixed" two-point vertex func- 
tion r^,, as depicted in Fig. ^| At the normalization 
point q = lj = 0, r=f,wc find 



1 tpij> — ~Mo 



I - 



uo(s + s' ) 



AD 2 
8D* 



f 



(k 2 

1 



2\2 



p 2 ) 



2s a s' + u Q (s' 



p 

so) 



p 2 f 
1 



P 



(3.36) 



Notice that the integral in the last term of ( 3.36j ) diverges 
in d = 2. In the same way as in the shift of the percolation 
threshold in DP, see Eq. ( 3.12 ), we may take care of this 
UV divergence by means of an additive renormalization, 
and then multiplicatively renormalize the UV poles in 
d = 4. Thus, defining the dimensionless renormalized 
coupling 



(3.37) 



with the associated flow function = Kd K ln( fj, / (mq), we 
have 



Moc — 



2s s + u (s + so) 



4A, 



f 



p K -(- p 



(3.38) 



and 



[ZqpZ^ 



) 1/2 Z, = 1 



u (s 



4D 2 
8Dg 



A k2 
I 



p 2 ) 2 



p (k 2 +p 2 ) 



2\3 



(3.39) 



We will lat er se e that the prefactor multiplying the shift 
/ioc in Eq. (3.38), which involves the various mixed three- 
point couplings, actually vanishes in an appropriate pa- 
rameter subspace containing both emerging fixed lines, 
see below. In principle, additional additive rcnormaliza- 
tions would be required to render t^rS and d u T^K 
UV- fini te [19] . These counterterms, to be added to the 
action (p72), would be of the form 



J d d x J dt(p(Ad t 



(3.40) 



However, as both A and B are again proportional to the 
prefactor of the integral in Eq. (3.38), they all vanish 
at the fixed lines to be discussed later. A subtle point 
which can be raised concerning these counterterms is the 
stability of the scaling behavior of the theory (in other 
words the stability of the fixed lines to be der ived later 
on) against the intr odu ction of a term like Eq. (3.40) into 
the original action (3.2). After this paper was submitted 
for publication, Jansscn |2S| has shown that indeed the 
scaling behavior is unaffected by the introduction of such 
terms, and thus /iq is the only mixed coupling constant 
that needs to be introduced. 

Note also that the fi nal di agram in Fig. [| [see the sec- 
ond lines of Eqs. ( ^36| ) and ( jo^ )] is UV-finite in d = 4, 



and so for the moment we shall neglect it in a minimal 
subtraction scheme (this is, however, a somewhat sub- 
tle point in the act ive ph ase, which will be discussed in 
more detail in Sec. Ill F ) . Certainly, for /j, <C 1, i.e., in 
an additional expansion in the transmutation rate /z, this 
diagram is suppressed as compared to the other contri- 
butions. We then find 



1 



M (SO 



8^ 



AD 2 



A 



(3.41) 



and after defining g — s/D and g' — s' /D, we have 

C M = -2 - 2v + V^(g + g') 



(3.42) 



where in the second line we have inserted the DP fixed 



point ( |3.35[) 
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FIG. 5. The "mixed" three-point vertex functions (a) 
Vtp^p, and (b) T^^,, to one-loop order. 



FIG. 6. The "mixed" three-point vertex functions (a) 
T^iiii, and (b) F^^^ to one-loop order. 



An inspection of Eq. (3.36) now shows that, in order to 
compute the renormalization of /xo, we must first consider 
the renormalization of the various mixed three-point cou- 
plings, 

s = ZsSoA 1 / 2 ^ 2 , s' = Zjs'oAWk-/* , 

a = Zfs A d /2 ^' 2 , s' = Z-^A^k-^ . (3.43) 

Evaluating first the renormalization of s', which can be 
calculated from the diagrams shown in Fig. ^](a), we find 



- <P<P>Pr ~ „l/2 „ „ 



U S S 



-u s + 



+ UqSq + U 

1 



s' J 2D 2 J p ( K 2 +p 2 ) 2 



(3.44) 



Inser ting the appropriate expressions for Z^ and Z v from 
Eq. ( |3.8| ) and using J p (n 2 + p 2 )~ 2 — A d n~ e /e, we end up 
with 



Z s . = 1 



bu 2 


UqSq 


u s a s 




l&D 2 


2D 2 


9 s' D 2 




u s' Q 


u Q s s' Q 


s 2 ;' " 


A d K~ £ 


2D 2 


9s' D 2 


9 s' D 2 


e 



(3.45) 



Similarly, with the aid of the diagrams shown in Fig. ||(b) 
and Fig. 6, we can compute the Z factors for the other 
mixed three-point couplings, 



Z s = l 



16D 2 



UqSq 




AD 2 


4s D 


u s' Q 


A d K~ e 




e 



(3.46) 



Z s = 1 



5^o 
WD 2 



UqSq 

2D 2 



2D 2 



mq£q£o 
2S D 2 

= 1 + 



UqSqSq 

) 2 



;/2„ 
s s 



2S D 2 2S D 2 



16D 2 
u s' 



AD 2 



4s' D 2 
A dK - e 



A d K 



AD 2 



(3.47) 



(3.48) 



We note in passing that, in principle, a product of the 
quartic vert ices [which we previously discarded from the 
action (3.1)1 and fiQ might also enter the renormaliza- 
tions of the three-point functions. However, we have 
checked that these additional couplings all have negative 
RG eigenvalues and are therefore irrelevant. 

With the definitions g = s/D, g' = s'/D, g = s/D 
and g' — s'/D, it is now straightforward to compute the 
RG P functions for these new variables. Since we have 
P a = nd K g = g(Cs - Cd) etc., where ( s = K<9 K ln(s/s ) 
etc., we find to one-loop order 



A? = + \/gS ( + -9 



e ( 1 



(3.49) 



v = Jhtf + d + J^/ig + g') + \g 2 ~9' , (3-50) 



Pa = \h9(~9 + ff) + \hsfiff +9) + t;~9 U 9 > (3-51) 



(3.52) 



where we have already inserted the one-loop fixed point 
value for v* = [(u/AD)*] 2 . Notice the symmetry of the 
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above RG j3 functions with respect to exchanging g <-> g' 
and g' <-> g. We now search for fixed-point solutions of 
the above equations where (3* — [3*, — /3| = flp =0. 
Using Eqs. (3-4E) and ( 3.52| ) to eliminate g and g' from 
the remaining two f3 functions, the above system of equa- 
tions can be solved exactly. After some tedious algebra, 
we find two fixed-line solutions, the first one being 



9 



g , 9-9 



2.g* , 



g* 2 = 2J e -(g*+g>*) 



(3.53) 



Computing the eigenvalues of the stability matrix 



<9/3 7 
85 



with M,{S} = {g,9',9,g'} 



(3.54) 



yields the eigenvalues 0, 0, — e/3, — e/3. Hence, this first 
fixed line, which includes the Gaussian fixed points for 
the mixed three-point couplings, is unstable. The second 
fixed line is given by 



(3.55) 



9 , g +g = 2W- 



* 2 = 2W^(s*+s'*) 



with stability matrix eigenvalues 0, e/3, e/3, 4e/3. Hence 
this second fixed line is stable. Notice also that both of 
the above fixed lines satisfy the condition 



-g g 



(3.56) 



which ensures the cancellation of the strongly singular 
(UV-divergent in d = 2) diagrams for the renormalization 
of in Eq. ( p6] ). 

We can n ow in sert t he values of g* +5'* at th e two fixed 
lines ( |3.53 ) and Q3.55 ), respectively, into Eq. ( 3.42| ), and 
thus obtain 



C 



-2 + e/6 + 0(e 2 ) (stable line) , 
-2-e/6 + 0(e 2 ) (unstable line) 



(3.57) 



We are now in a position to derive the scaling form (2.3C) 
for the B species density, postulated in the previous sec- 
tion. We begin by writing down the renormalization 
group equation for the renormalizcd field (<p_r), where 
the angular brackets denote averaging with respect to 
the full action S mc 



d „ d d d 

ok or oD ov 

• <*°Tg • ^ ' -4 • (3 ' 58) 
d 1 \ 

and where we have used the notation {g} = {g,g',g,g'}. 
Defining the dimensionless field <p as 



(<Pr(k, t, D, v, {g}, fi, x, t)) = n d/2 

xip(T,v,{g},fi/D,Kx,K 2 Dt) , (3.59) 



the solution of Eq. (3.58) when the couplings v, {g} have 
run to their fixed-point/linc values is 

r, D, {g}, n, x, t)) = n d ^ d -0/z 
x$(riG,v*,{g*},(ji/D)ft-&, 

kx£,k 2 DU 2+Cd ) . (3.60) 



Inserting the matching condition I = |r| 1 ' ^ , and drop- 
ping the i>, {g} couplings, we obtain 



(<p R {K,T,D,H,X,t)) OC \ T \-{d-C;)/2C 



(3.61) 



-«:-c D )/c 



K 2 Dt\T\-v+&yc 



-i/C 



Identifying ft = = -{d-C 9 )IK*, "X = 1 
— (2 + Cz>)/ C* an( i defining the crossover exponent as 4 
(Q - C D )/C, we have 



(3.62) 



in agreement with the scaling hypothesis ( [2.30| ) postu- 
lated in the previous section. Using our earlier results 
for the C* functions, we find 



_ / 1 + <3(e 2 ) (stable line) , 

1 + e/6 + (9(e 2 ) (unstable line) 



(3.63) 



Notice the absence of O(e) contributions to <j> at the sta- 
ble fixed line, which is due to remarkable cancellations. 
This of course also implies that there are no logarithmic 
corrections to the crossover exponent <j) in d c = 4 dimen- 
sions. 

The final step in this calculation is now to compute 
the exponent ft- Unfortunately, in order to do this, we 
must first understand the behavior of the scaling function 
( |3.62 ) in the active phase. As we shall see, it contributes 
non-trivial corrections to the exponent ft at O(e). Hence 
it is not sufficient simply to match the scaling function 
to that calculated in mean-field theory |^3| — such a 
procedure would miss these O(e) corrections. However, 
before dealing further with this active phase calculation 
for the B species, we first discuss the simpler problem of 
pure DP in the active phase, which also applies to the 
first level (the A particles) of the coupled DP problem. 



D. DP field theory: Active phase calculation 

In this section, we will review the one-loop calculation 
of the expectation value of the field in the active phase 
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for the case of a single field, i.e., the case of (decoupled) 
directed percolation. In this way we can obtain an ex- 
pression for the critical exponent f3\ = j3 . 

W e sta rt with the action for a single field (see 
Sec. |fflAl) 



Si 



DP 



d d x dt{4> [d t + D (r -V 2 )] Vo- 



. (3.64) 



In the active phase (ro < 0) the expectation value of the 
field V'o is non-zero, and we define a shifted field V> c o by 

ipo = vo + V>cO , (3.65) 
to obtain the new action 

S'dp = d d x dt {Mdtipco - D V 2 ip c0 + D roipda)+ 



(3.66) 



1/2 u v 




FIG. 7. One-loop diagram for {ip c o) in the pure DP active 
phase calculation. 



We see that all the poles in e cancel out as they should. 
At the fixed point (u/D)* = 2(e/3) 1/2 , we finally find to 
leading order in e, 



(il>n) = 2 ( - ) A l J 2 K d /\l + e/6) r (l - ~ lnr 



(3.72) 



We now fix vq by equating the coefficient of tpo to zero. 
Thus 



t'o 



2D r 2D \r \ 



Uq 



«o 



(3.67) 



where vq is the classical (mean-field) value of the expec- 
tatio n val ue of the field. Substituting these values into 
Eq. (3.66), we obtain 



Upon exponentiating the logarithm, assuming that this 
can be done unambiguously, we see that 



(iPr) = n A ~ r 1 -^ 6 



(3.73) 



which i mplie s (3\ = /3 = 1 — e/6 + 0(e 2 ), as already cited 
in Eq. ( |3.22 ) , where the scaling relation /3 = ux (d + z — 
2 + rj±_)/2 was employed. 



St 



DP 



d d x / dt [Mdt^co - D V 2 ^ c0 + D \r \ipco)+ 



(3.68) 



In the following we will define tq = —tq. From Eq. (3.65) 
it follows that 



(ipo) = v + (ipco) 



(3.69) 



There is only one diagram contributing to the expecta- 
tion value {ipco) to one-loop order, and it is depicted in 
Fig. @. We thus find, using dimensional regularization, 



1 



2Dqt _ 1 «o_ 

u 2 D J p p 2 + t 
2D t u T(l - d/2) 



uq Do (8ir) d / 2 



(2r 



.d/2-1 



(3.70) 



Using the relations betw een t he ba re and renormalized 
quantities given in Eqs. (|3 . 7|) — ( 3.11 ), the last expression 
yields 



( 1 + 



^ D 

u 2 

AeD 2 1 - e/2 



(3.71) 



E. Coupled DP field theory: Active phase 

We now proceed to discuss the active phase calculation 
for two coupled fields, as represente d by t he action S mc — 
S e g + AS given in Eqs. (§J) and ( |3~27| ). We will again 
restrict our discussion to the case where ta = rs = to- 
Based on the mean-field analysis of Sec. 0, we anticipate 
that the expectation values of the fields are different from 
zero when tq < 0. Hence we define the shifted fields 



IpQ = VA0 + IpcO 
1p0 = tpcQ , 

(pa = vbo + <Pc0 

<P0 = <fc0 ■ 



(3.74) 
(3.75) 
(3.76) 
(3.77) 



The new action expressed in terms of the shifted fields 
will now involve terms linear in ip c o and (p c o- Equating 
the coefficients of these terms to zero fixes the constants 
vao and vbo to be the classical (mean-field) values. These 
are determined by the equations 



1 2 

D r VAo + 2 U o v ao = i 

n l 1 2 

D r VBo + 2 W o«bo ~ HoVao 



- t s v A o 



(3.78) 
s' vaoVbo = . 
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FIG. 8. New two-point vertices in the active phase coupled 
DP calculation. 



The solutions of these equations are 
2£>oM 



VAO 



VBO 



»0 



/4/4 A)|n)| 



O(|ro|) 



(3.79) 
(3.80) 



1/2 u v A0 




(e) 



1/2s' v A0 



(f) 



s V A0 




(g) 



1/2u v A0 





(h) 



The fields now have new masses given by 

D rAo = D a r a + u vao = D \r \ , (3.81) 

D r B o = D r Q + u v B q + §' Q v A0 (3.82) 

~ ^A^D \r \ + O(\r \) . (3.83) 

There are also various new two-point vertices as depicted 

in Fig. ||, whose corresponding terms in the action are 
given by 



d d 2 



dt 



l;u VAoiplo - t;{u vbo + s' a v A0 )(pl Q - 



-(M0 - 8qV A - SqVbO^cOiPcO - so^ao^'coV'co] ■ (3.84) 

We can now proceed with the calculation of th e exp ecta- 
tion value of the second field 1^0, by using Eq. ( 3.76 ) and 
computing (ip c o) to one-loop order. The full calculation is 
unnecessarily complicated, and universality dictates that 
we can calculate the critical behavior for any point on 
the stable fixed line. However, we note that the unstable 
and stable fixed lines do not necessarily have to give the 
same critical behavior, and indeed in our case we will find 
that they do not. 

For the stable fixed line, we have actually performed 
the calculation in two different ways. Both methods yield 
identical results, and this provides us with an important 
extra check on our methods. In the first approach we 
have chosen a subspace of initial parameters 



•so 



5 0' °0 



•so 



(3.85) 



with the stable fixed point (u/D)* = 2y/e/3, g* 
—2g'* = —2g* = \/e/3. In the second approach we have 



chosen 



Sn = s = s' = , 



(3.86) 



FIG. 9. One-loop diagrams for the computation of (<j> c o) 
in the coupled DP active phase calculation. 



and with the stable fixed point (u/D)* = 2y / e/3, g* = 

(s/D)* = 2^/7/3, 9'* = 9* = 9'* = 0. 

No te th at bot h th e parameter subspaces given by 
Eqs. |[8|) and ( ^86| ) are closed under renormalization 
flows, as is evident from a clo se insp ec tion of the RG ( 3 
functions given by Eqs. (§T9|), ( ^50| ), |[|l]) and ( ^52| ). 
The first choice is the more natural one for investigating 
the multicritical point, since when Eq. (3.85) is satisfied 
the original action S mc = S c s + AS is invariant under a 
generalization of the usual DP 'rapidity reversal', namely 



-tp (x,-t) , 
-<Po(x,-t) , 



(3.87) 
(3.88) 



and under renormalization this symmetry is pr eserv ed. 
However, the calculation using the condition ( |3.86 ) is 
somewhat simpler. In both approaches the unstable 
(non-trivial) fixed point is given by (u/D)* = 2-\/e/3, 
and g* = g'* = g* = <?'* = 0. 

The diagrams contributing to the expectation value of 
((fico) to one-loop order are depicted in Fig. |[ We will 
outline first the calculation using the second approach 
described above. In that case diagrams (e) through (h) 
vanish and one has to consider only the first four dia- 
grams. We th us find to one-loop order, after implement- 
ing Eq. (pc|), 



(ipo) = vbo - h(r B o) - 



uIvao^I 



2D 



r B0 



2D r B o 
h(rA0, r B o) 



2D^r AQ r B o 



h(rAo) 



D r B o 



h(rBo\ r A o,r B o) 
+ ••• . (3.89) 
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With the use of dimension al re gularization, the various 
integrals appearing in Eq. ( 3.89 ) are given by the follow- 
ing expressions 



h(a) 



h(a,b) 



1 



1 



A d 



2D J pP 2 + a D e(2 - e) 

i r i 



A d 1 



a 1 " 6 / 2 , (3.90) 
(3.91) 



D*e(2-e) (a -6)= 



a !_ 6/2 , & l_ e /2 



l-e/2 



h(b;a,b) 



±DlJ p (p2 +6)(p2+ a±6) 



-4, 



2Z? 2 e(2 - e) 



a - 



a - 



a - b V 2 
2& 



e/2 



6-e/2 



6 — a 



We are interested in obtaining the e expansions of the 
above expressions in the regime where b 3> a (since 
tbo tao as rg — ► 0). In that case we find 



A d K 



2D e 
A d K~ e In 2 



l--ln4 

2 K 2 



(3.93) 



H In 2 1 In 



7i(a) 
/ 2 (a,b) 
I 3 (b;a,b) 



Let us postpone for the time being a full discussion of 
the term containing I2, i.e. the contribution from dia- 
gram (c) in Fig. ^. Notice that this term is ultraviolet 
finite, i.e. it docs not i nvolve a pole in e. However, after 
using expression ( 3.93 ) it appears to contribute a (nega- 
tive) constant to the expectation value of (p. This seems 
problematic and a full discussion of this point will be 
given in the next section. 
We now define again 



to = -(ro - r 0c ) , 



(3.94) 



where to first order in standard perturbation theory we 
have 



roc 



1 



Jp V 



(3.95) 



This integral vanishes in dimensional regularization, and 
hence in the following calculation we will ignore tq c , 
as was done in the pure DP active phase computation. 
When using the relations b etw een th e bar e and renor- 
malized quantities in Eqs. (3.7) and ( 3.37 ), we find to 
leading order in an expansion in t: 



gyjlDr x/2 d /2 7 i/2 7 7-1/27-1/27-1/2 



4D 2 , 



■(1 + 6/2) 1 



iln^) 
4 D I 



u 



+ F7^r;( 1 + e / 2 ) 1 



8L> 2 e 
us 



lnr 



(3.96) 



sl)2 e (l + ^-eln2/2)(l-£ln^) 



A check reveals that all the poles in e cancel, as they 
should. At the fixed point given by (u/D)* = (s/D)* = 
2i/e/3, we find 



(3.92) + J) 



6J V D 

I 1 - -^H^r/D) - ^7 lnr 



(3.97) 
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Assuming that only one such scaling term becomes gen- 
erated, we may now exponentiate the logarithms, and see 
that 



l/2-e/24 



r W6 ( r 



from which we infer 



0=1 + O(e 2 ) , 



(3.98) 

(3.99) 
(3.100) 



where <f> is the crossover exponent, a result consistent 
with that derived previously in the inactive phase. We 
also see that in the definition 



(<Pr)=t*-<p(t-*h/D) 



(3.101) 



the scaling function behaves for large argument as 



_l/2-e/24 



(3.102) 



a result that differs from the mean-field behavior at O(e). 
At the up per critical dimension d c = 4, a comparison 
with Eq. ( [3.26 ) for the first hierarchy level suggests the 
logarithmic correction 



(^-rV^-lnr) 1 ^ 



(3.103) 



Let us also discuss briefl y the calculation using the 
approach given by Eq. (3.85) above. In that case we have 
an additional four diagrams (e)-(h) to consider. However, 
we will see that diagrams (f)-(h) contribute to higher- 
order terms in an expansion in r, beyond the leading 
behavior described above. Diagram (e) gives an equal 
contribution to diagram (d), but since the value of s at 
the corresponding fixed point is half of what it was in the 
previous calculation, the final result is exactly the same. 
To verify these claims we observe that diagrams (e) to (h) 
cont ribut e the following additional terms to the r.h.s. of 
Eq. (p9b , 
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uaSoVAofJ-o 



2D r B0 



h{rAo;rBo,rAo) 



uqs vaq 
2D r B0 



h(r B o) 



S^VAO 



h(r A o,r B o) 



UqSqVAO 



D r B Q 2D r B0 
The only new integral I 4 is given by 

h{a,b) = h((a + b)/2) 



h(rAo) 



(3.104) 



(3.105) 



It is now easy to verify that the last three terms give a 
contribution of O(ro) and higher which can be neglected 
in comparison with the 0(^/tq) terms which we have 
kept. The I3 contribution in the first term is given by 



I 3 {a;b, a) = 



2D 2 e(2 - e) 



b + a ( a + b 



2a 



-e/2 



-e/2 



4Dge 



1 



In 2 



- b 



1 - -ln6 
2 



(3.106) 



and thus diagram (e) yields essentially the same contri- 
bution as diagram (d). Since the value of (s/D)* is now 
one-half of its value in the previous approach, our final 
result follows. 

Finally, let us discuss the behavior at the unst able, 
non-trivial fixed point given by (u/D)* = 2 \/ e/ 3, and 
g* = g'* = g* = g'* = 0. Returning to Eq. ( pl96| ), with 
the last term set equal to zero, we have 

(<PR) « J% (l - ^ H^/D) - ± lnr + . . 

(3.107) 

to leading order. Exponentiating the logarithms (again 
with the assumption mentioned above) we find 



\/2-e/\2 



(3.108) 



(3.109) 



and thus in this case fa = l/2-e/6+0(e 2 ), <j> = l + e/6+ 
0(e 2 ) and the scaling function p(x) ~ j.i/2— e/i2+0(e ) f or 
large values of the argument. 



F. Technical difficulties 

We now return to diagram (c) of Fig. 9. Substitut- 
ing the expression for l2(a,b) fr om E q. (3.93) into the 
corresponding expression in Eq. ( |3.89| ), we find that the 
contribution of diagram (c) to the expectation value of 
(<Po) is: 



ulA d n e V A Q o 



In 2 /! 
T6~D' 

(3.110) 



to leading order in e and in an expansion in powers of 
yfr. This is in contrast to our expectation (backed up by 
our simulation results in Sec. IV) that (ip) should vanish 
at the transition as t — > 0. Notice that this diagram 
is ultraviolet-finite and hence there are no poles in e. 
On the other hand, the loop integral 12(0,, b) is infrared- 
divergent for any d < 6, which leads to the behavior 
l2(a, b) cx 1/b for d = 4 when a — > 0. One can argue that 
for d > 4 the prefactor of this diagram will vanish due to 
the fact that u flows to the Gaussian fixed point u* = 0. 
But for d < 4 we seem to have a problem, the origin of 
which is obviously the appearance of the strongly relevant 
parameter (i as an effective coupling (two-point vertex) 
in the perturbation expansion. 

This difficulty is somewhat reminiscent of the random- 
field problem, where infrared-divergent diagrams in per- 
turbation theory lead to a shift in the upper critical di- 
mension in a tp A theory from 4 to 6 | p5[ . This is due 
to the extra propagators in a given loop as compared 
to the pure case. One might perhaps argue that be- 
cause of the unidirectionality of the interaction between 
the two fields ip and ip, the ip field acts as a spatially 
and temporally correlated random field from the point of 
view of the ip field, since there is no backwards feedback, 
meaning that the <p field has no influence on the i/j field. 
However, in contrast to the random field case, where the 
random field is taken to be of constant variance as one 
approaches the transition temperature, the expectation 
value of the ip field vanishes as one approaches the mul- 
ticritical point, and this seems to soften the effect of the 
infrared problem. Another marked difference is that "0 is 
not a quenched random variable in the traditional sense, 
since it is displaying strong temporal fluctuations in the 
multicritical regime. Since the above infrared-divergent 
integral becomes tamed at d = 6, one might think that, 
similar to the random-field problem, one might control 
its divergence by a dimensionality shift. However, we 
have not been able to construct a sensible field theory 
for this problem by choosing d — 6 as the upper critical 
dimension. Actually other infrared-divergent diagrams 
can be found at higher loop orders, which are even more 
divergent than the diagram just mentioned, as in, e.g., 
Fig. [lCj(a). IR-divergent diagrams also appear in the ex- 
pansion of other vertex functions, like the pip 2 vertex 
to one-loop order, where one can easily construct dia- 
grams with one or two insertions of the ip field into the 
(^-triangular diagram, see Fig. 10(b),(c). 

Another place were a similar diagram appears is in 
the "mixed" vertex function T^,, see the last diagram of 



Fig. 4 and Eq. ( 3.36 ). If we evaluate this vertex function 
at a ge neral "temperature" r, rather than at r = 1 as 
in Eq. ( |3.3fi| ), we find that its contribution in d = 4 di- 
verges like 1 / t as t — > . From the active phase side, the 
same diagram diverges like 1 / ^JJIf. This prevents one 
from defining the renormalized \i parameter at critical- 
ity as the value of this function for t = 0, at least not at 
q = to = 0. One can attempt to absorb this divergence by 
a finite (non-ultraviolet divergent) renormalization addi- 
tion to Z^, which will then become temperature depen- 
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FIG. 10. (a) An IR-divergent diagrams at higher loop 
order, contributing to the expectation value of (p. (b) 
IR-divergent diagram contributing to the three-point vertex 
function with one insertion of tp 2 . (c) IR-divergent diagram 
contributing to the three-point vertex function with two in- 
sertions of tp 2 . 



Thus one would define such that the renor- 



dent. 

malized ITw, evaluated at q = lu = would be finite as 
r — * 0. If this is done in the active phase we have verified 
that it cancels the contribution of the problematic con- 
tribution to (tpo) against the corresponding term in Z^. 
However, this procedure appears somewhat artificial and 
not entirely satisfactory since it requires a temperature- 
dependent renormalization constant. Also, to establish 
such a renormalization scheme, one would have to check 
that this procedure works properly to higher loop orders 
as well, where even more IR-divergent contributions are 
present. 

In the absence of a truly satisfactory resolution of the 
infrared difficulty, one might argue that its effect is to 
render the e expansion derived in the previous sections in- 
valid. Yet, this conclusion is too far reaching, since even 
if a well-defined field theory that describes the asymp- 
totic critical behavior of all hierarchy levels is problem- 
atic (and may even not exist), our results are still likely 
to be valid for a range of the parameter r close to the 
transition point, as will become clear from the following 
reasoning. 

We notice that the problematic IR-divergent diagrams 
are proportional to higher powers of the inter-species cou- 
pling \i. For example, the IR-problematic diagram men- 
tioned at the beginning of this section gives an overall 
contribution proportional to fj,, whereas the other dia- 
grams contributing to (tp) to one loop order (which were 
mentioned at the end of the previous section) are pro- 
portional to y/JI or ^JUa/j,. For small values of \i these 
problematic diagrams are therefore suppressed. Now \x 
is a relevant coupling in the RG sense and thus the ef- 
fective (running) coupling increases as one goes deeper 
into the critical region. But by starting with an ini- 
tially small value of /i, one can increase the size of the 
region in which ji remains small. In this intermediate re- 
gion the IR-divergent diagrams can still be neglected and 



the scaling results obtained in the previous sections are 
valid. Thus, our theory predicts that for small /x, as r 
is decreased, one can observe a scaling regime where the 
critical behavior is characterized by the universal expo- 
nents calculated in the previous sections. In particular 
the exponent /?2 should be observable. Ultimately, as 
t becomes very small the field theory may break down. 
However, we cannot exclude the emergence of another 
non-trivial asymptotic scaling regime deep in the criti- 
cal region. The ensuing scaling behavior might possibly 
be extracted by isolating the structure of the leading IR 
divergences to all orders and then by a resummation of 
the ill-defined perturbative expansion in ji. However al- 
though this possibility exists, it cannot be substantiated 
at this point by more rigorous arguments. What is more 
likely to happen is that a non-universal crossover behav- 
ior ensues (non-universal as a consequence of the break- 
down of the e expansion and hence the RG construction) , 
which eventually terminates in asymptotically decoupled 
DP behavior. In other words, the B species are no longer 
slaved to the A particles, but rather behave indepen- 
dently such that their density vanishes with a power f3±. 
This is supported by the simulation results, as will be dis- 
cussed in a later section. Probably, at sufficiently large 
times, the discrete and finite number of A particles al- 
ready vanishes in the interior of comparatively large do- 
mains. Consequently the B particles, which were actu- 
ally generated previously through the reaction A — ► B, 
interact and annihilate as if they were independent and 
decoupled from the A species. Further discussions of this 
issue in the light of the simulation results will follow at 
the end of the paper. 



G. Diagonalized theory 

We now return to another of the appro aches to the 
coupled DP problem mentio ned in Sec. Ill A, namely that 
of diagonalizing the action (3.2) to remove the quadratic 
cross-term linking the tpo and <po fields. If we apply the 
transformations 



*o = <A) + 
#o = tP Q - 



Mo 



D (r A - r B ) 
Mo 

D (r A - r B ) 



ipo , 



(3.111) 



then the action (3.5) is transformed to 
Sdiag = Jd d x J dt {*o [d t + D (r A - V 2 )] Vo- 



.|[^ -#oVo 2 ] 
u 



(3.112) 



-cp [d t + D (r B - V 2 )] $ - y [^o*o ~ ^o$o] 

S S' 1 



Here we have defined 
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(3.113) 



Notice that this new action ( fc.llS ) has precisely the form 
of our original full action S mci except that the quadratic 
cross-term which linked the ipo and ifo fields has been 
eliminated. This renders the diagonalized model a spe- 
cial case of the very general quadratically coupled DP 
processes studied by Janssen |l5). In addition, we re- 
mark that this diagonalizing transformation breaks down 
when r A — r B , consistent with our earlier identification 
of novel multicritical behavior for r A = r B = r — > 0. 

since the action ( p. 112; ) is very similar to S mc , we 
can quickly determine its ensuing scaling behavior. The 
mixed three-point vertices hav e exa ctly the same fixed- 
line structure as derived in Sec. Ill C . Thus the diagrams 
which would have generated a quadratic cross-term under 
renormalization (i.e., the d — 2 UV-divergent diagrams 
in Fig. |J) cancel out at both of the fixed lines (at least to 
one-loop order). Hence, to this order, the DP parts of the 
action for the ip and $ fields are entirely unaffected by 
the presence of the mixe d thre e-point vertices, generated 
by the transformation ( 3.111 ). In that case we expect 
pure DP behavior on the transition lines away from the 
multicritical point, as we had earlier anticipated, and as 
shown on very general grounds in Ref. VM. 



H. Crossover theory 

Outside the multicritical regime, the mean-field ap- 
proximation suggested that we should expect ordinary 
DP critical behavior for the B species, if either t a > 
and t b — > 0, or t b > and t a — > 0, see Fig. pi If one 
starts near the multicritical point where both t a — > 
and t b — > 0, then at first the multicritical density ex- 
ponent /?2 should be observed, eventually crossing over 
to the ordinary DP exponent (3\, and similarly for the 
exponents cti, etc. In the first crossover region t a > 
and t b — > 0, indicated by the dashed arrow B in Fig. g, 
this scenario is rather obvious on the mean-field level, 
valid for d > d c : Here, the density n A = 0, and there- 
fore the coupling between the different hierarchy levels 
vanishes. Our stud y of the diagonalized theory in the 
preceding Sec. [II G established the DP character of the 
ensuing active/ absorbing transition at t b = even below 
d Cl when fluctuations become dominant. 

The more interesting case is the situation for t b > 
and t a — > 0, where there exists a DP phase transition 
for the A particles. For [i AB > 0, we argued above that 
the B species becomes "slaved" to the critical A species, 
and is driven to a non-equilibrium phase transition itself 



as the A control parameter t a — » 0. In order to describe 
the ensuing crossover scenario, we apply the generalized 
minimal subtraction scheme of Ref. (26) , where we retain 
the parameter t b in the RG flow equations, and use t a = 
1 as the normalization point. The Z factors then become 
functions of the non- linear couplings and t b . In order to 
simplify the calculation, we once again use the reduced 
parameter space where g' = g = g' = 0. However, now 
we must distinguish between the non-linear DP coupling 
v for the A species, and the one for the B particles, which 
we denote with v' . The remaining important Wilson RG 
functions then become 



C TA = -2 + 3v , 



3v' 



(i + r B y+^ ' 



Cu - Cd = -2 - v 



v'g 



{\ + r B y+^ 1 



(3.114) 
(3.115) 

(3.116) 



and the RG (3 functions for the couplings v, v', and g 
read 



v =v (-e + 12v) 
f3 v > = v' \-e 



12v' 



Pg=9 \ + 



(l+T B )!+ £ /2 

v'g 



(1 + t b ) 1+£ / 2 



(3.117) 
(3.118) 

(3.119) 



Of course, the A species DP fixed point remains v* = 
e/12 with C* A = -2 + e/4. However, the effective B 
coupling v' (l)/{l+T B (£)] 1+e / 2 — > asymptotically as £ — > 
(since v'jtj ~ £~ e ) and h ence t b (£) ~ £~ 2 according 
to Eqs. (|3.118| ) and (|3.115| ). Thus /3 g -> -(e/3)g, and 
g(£) -> as well. Consequently, (* - (* D = -2 - e/12, 
and we recover 



= (C-G)/C-i + e/6 + o( e 2 ) 



(3.120) 



i.e., the result given in Eq. (3.63) at the unstable fixed 
line. Notice that here the critical exponents for the B 
species are given by v±_ = —1/ C* e tc., and thus take on 
the usual DP values fl3.19]) to (§J2§). 

In the active pha se, we ma y con firm this by direct cal- 
culation, as in Sees. [II D| and III E. From the diagrams in 
Fig. |9| for finite t b , and in essentially the reduced param- 
eter space, only graph (b) survives. After inserting the 
p c shift, and multiplying by the appropriate renormal- 
ization constants (in the generalized minimal subtraction 
scheme), the equation of state in terms of the renormal- 
ized quantities assumes the form 



(<Pr)tb\ta\ 



2v^^(^.) 2 



In I ta 

2 2 1 



(3.121) 



At the DP fixed point v* — e/12, we exponentiate the 
logarithm on the left-hand side to Ir^l 1 " 
use (ipR) ~ t^ 1 , which leads us to 



e/6 



r^ 1 , and 
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(3.122) 



which obviously generalizes the mean- field result j2~2^ ) 



I. Higher hierarchy levels 

We end this section on applications of the renormal- 
ization group to coupled DP with some remarks on the 
behavior of higher hierarchy levels i > 2. As pointed 
out already in Sec. ||, while generically the transitions 
from the active to the absorbing inactive phase of particle 
species i are of the DP universality class, with the critical 
density exponent (3\ , one hnds the two-level hierarchy ex- 
ponent fa whenever the first two critical points coincide, 
ta =tb> However, further special multicritical behavior 
appears at higher levels if additional become equal. In 
mean- field theory, one has ft = l/2 l_1 for r\ = . . . = r.i. 
In the two-level hierarchy (i = 2) multicritical regime, 
one finds a strong downward re norm alizati on of ft and 
ft due to fluctuations, see Eqs. ( |3.22| ) and ( 3.99 ). 

In order to see what happens at higher hierarchy lev- 
els, let us briefly consider the three-species coupled DP 
process in the vicinity of the multicritical regime, and 
for small transmutation rates. Notice that in this situa- 
tion fluctuations w ill no t only generate the vertices cor- 
responding to Eq. ( 3.27 ) on each adjacent level, but also 
the indirect transmutation A — > C, as well as all possible 
three-point vertices unidirectionally coupling the levels 
(and in principle additional higher order non-linearities, 
which however turn out to be irrelevant). In order to 
simplify the analysis, we merely use the reduced parame- 
ter space analogous to keeping only the new vertex sq in 
the two-level process. This leaves us with the following 
effective action (with the fields (j>o and </>o describing the 
particle species C), 



SeS z 



d d x / dt 



{-00 [d t + D (r A - V 2 )] </> 



-y (V>oVo - ipQipo) + 

1t 



(3.123) 



+ip [d t + D (r B - V 2 )] <po - — (^o^o - VWo) ~ 
+4 [dt + D (r c - V 2 )] 0o - y 1 (4>Uo - Ml) - 

-/i' 4> ip - t 4> ip ip Q - 

Clearly, if we just consider the B/ C rea ctions, then a 
fully analogous calculation as in Sec. Ill C yields 



Cm' - Cd = -2 - v 



rvh 



(3.124) 



where h = t/D. Hence with v* = e/12, the crossover 
exponent at the associated two-level multicritical point 
is 



^ = (C -&)/£ = 1 + 0(0, 



(3.125) 



computed at the stable fixed point h* = 2y / e/3 (which 
is the remnant of the stable fixed line in the reduced 
parameter space). We may wonder now if the unidirec- 
tional, sequential coupling of three hierarchy levels leads 
to a further novel crossover exponent associated with the 
transmutation A — > C. Thus, we compute the renormal- 
ization of the vertex function T^, which leads to 



Cn - Cd = -2 - v + V^(f + f) 



(3.126) 



where / = t/D and / = pfi'/Dp,. We now merely need 
the renormalizations of both / and /. In just the same 
manner as in the two-level calculation, one finds from the 
renormalization of T^^, 



(i f = /(-e/2 + 2v + yfif) = /(-c/3 + Ve/12f) , 

(3.127) 

which obviously has the IR-stable fixed point /* = 
2-v/e/3. The only really novel renormalization concerns 
the vertex function T^^, yielding 



Cp - Co = -e/2 -2v + ^(g + h + f) 



(3.128) 



After using — = V^{h — f — /), this leads to 



f3 f - = /[-e/2 - 2v + ^ (g + 2h - /)] = 



= f(e/3-^7/l2f) 



(3.129) 



Comparing with Eq. ( |3.127 ), we see that the non-trivial 
fixed point /* = 2\J e/3 is unstable, whereas /* = is 
stable for d < 4. Consequently, the three-species vertex 
T^^/j vanishes asymptotically, i.e., becomes irrelevant, 
which implies 



<^(C*-CD)/C = l + 0(e 2 ) 



(3.130) 



identical with the one-loop values of the crossover expo- 
nents <f> and <j)' . Coupling to an additional hierarchy level 
therefore does not introduce a novel crossover exponent 
in the multicritical regime, at least to O(e). We suspect 
that this is actually true to higher orders in e = 4 — d 
and for higher levels i > 3 as well. 

On the other hand, below the critical dimension d c — 4 
the density exponents fii are affected by the fluctua- 
tion corrections to the scaling functions, and ar e no t 
simply determined by a scaling relation like Eq. (2.37) 
p3| . In fact, were such a renormalization contribution 
from the scaling function absent, we would arrive at 
Pi = l/2 l ~ 1 — e/6 to one-loop order, with the O(e) correc- 
tion independent of the level index i. This would predict 
that near four dimensions, i.e., for any < e < 1, ft 
should bec ome n egative for sufficiently large i. The cor- 
rect result (3.99) for the second hierarchy level, however, 
shows that the O(e) correction is actually smaller than 
for the first level. Presumably, on each successive level 
the O(e) corrections are further reduced, such that all the 
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FIG. 11. Schematical illustration of a coupled hierarchy of 
three directed bond percolation processes in 1+1 dimensions. 
Activity percolates along the direction of time through bonds 
(simple arrows) which are open with probability p. The three 
subsystems are coupled by instantaneous transfer of activity 
(double arrows) to the corresponding site of the next subsys- 
tem. 



density exponents remain positive. A detailed computa- 
tion of /?3, which requires an explicit study of the active 
phase, already becomes a rather tedious affair, and we 
leave our discussion of higher hierarchy levels with this 
speculation. 



Coupled DP may be realized on a computer by si- 
multaneously evolving several directed bond percolation 
models which are coupled without feedback according to 
the principles outlined in the Introduction (see Fig. pL l[ ) . 
For simplicity, we work in the limit of infinite coupling 
strength fj, = oo, i.e., active sites in one of the subsystems 
instantaneously turn the corresponding sites of the next 
subsystem into the active state. Alternatively, we could 
also apply a finite coupling strength (probabilistic trans- 
fer of activity to the next subsystem), or explicit parti- 
cle transmutation A — > B,B — > C, . . .. Although these 
variants are expected to have the same critical proper- 
ties, their initial crossover times into the scaling regime 
are typically longer. Therefore we restrict the numerical 
analysis to the case of infinite coupling strength. 

In principle it would be possible to simulate arbitrar- 
ily many hierarchy levels. However, in order to reach 
the scaling regime, the particle densities have to be low 
enough. It turns out that already at the fourth level the 
particle density is rather high, which makes it extremely 
difficult to determine critical exponents. For this reason 
our numerical simulations are restricted to three hierar- 
chy levels. 



A. Numerical estimation of the critical exponents 



IV. NUMERICAL RESULTS 

In order to support our field-theoretical results, we 
study a unidirectionally coupled hierarchy of DP models 
using Monte-Carlo simulations. There is a large variety 
of DP lattice models that can be used for this purpose. 
One of the simplest and most efficient realizations is di- 
rected bond percolation on a tilted square lattice [^7| . In 
this model, neighboring sites are connected by directed 
bonds which are open with probability p and closed oth- 
erwise. Activity percolates through open bonds along a 
given direction which is usually interpreted as the direc- 
tion of time. Labelling different rows of sites by a discrete 
time variable t, directed bond percolation may be equiv- 
alently defined as a stochastic cellular automaton with 
parallel update rules mapping the system's configuration 
at time t probabilistically onto a set of new configura- 
tions at time In one spatial dimension, directed 
bond percolation is just a special case of the Domany- 
Kinzel cellular automaton [^8| which is known to be one 
of the most efficient realizations of DP on a computer. 
Another advantage of using directed bond percolation is 
the availability of very precise estimates for the critical 
percolation threshold p c (d) in d < 2 spatial dimensions. 
Currently the best estimates are p c (l) = 0.644 700 15(5) 
H for d = 1 and p c (2) = 0.287 338(6) @ for d = 2, re- 
spectively. To our knowledge, the percolation threshold 
for (3 + l)-dimensional directed bond percolation has not 
been estimated before. Using standard methods we find 
the value p c (3) = 0.13235(20). 



In order to estimate the critical exponents of coupled 
DP, we employ two standard numerical methods for sys- 
tems with phase transitions into absorbing states. On 
the one hand, we use steady-state simulations in the ac- 
tive phase in order to directly determine the exponents 
(3k- On the other hand, dynamical simulations |3l|] at 
the multicritical point render a set of dynamic exponents 
which in turn determine the exponents fj_ t k and v\\ : k- 

a. Steady state simulations in the active phase. On 
the multicritical line the stationary particle densities nk 
are expected to scale as nk ~ (p — p c Y h ■ By measuring 
nk in a sufficiently large system, it is therefore possible 
to directly estimate the exponents (3k- The accuracy of 
the results depends on the accessible range of Ap = p — 
p c in the simulation. In d = 1 spatial dimension the 
minimal value of p — p c is predominantly limited by the 
equilibration time T equ which has to be larger than the 
temporal correlation length ~ (p — p c )~ U] ^ k , whereas 
in d > 2 dimensions the main limitation is the system size 
]V max , which has to be larger than £j_ k ~ (p — p c )~ dw± - k 
(see Table |) . The measurements for d — 1 are shown in 
Fig. 12 1. From the slopes of the lines averaged over one 



decade we estimate the exponents /3j, whose values are 
listed in Table §. 

b. Dynamical simulations at the multicritical point. 
The most precise estimates for the critical exponents of 
systems with phase transitions into absorbing states are 
usually obtained by dynamical simulations plfl . Starting 
from an initial state with a single particle (active seed), 
the system evolves at the critical point and generates 
a spatio-temporal cluster of active sites whose size and 
lifetime is finite. Survival probability, mass, and mean 
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d = 1 


d = 2 


d = 3 


d = 4 - e 


Ap min 


0.0004 


0.0008 


0.0016 




■^V max 


4000 


100 2 


35 3 




Xequ 


10 s 


10 4 


10 3 




Pi 


0.280(5) 


0.57(2) 


0.80(4) 


l-e/6 + 0(e*) 


132 


0.132(15) 


0.32(3) 


0.40(3) 


l/2-e/8 + 0(e 2 ) 


Ps 


0.045(10) 


0.15(3) 


0.17(2) 


1/4 - O(e) 


Pop 


0.2765 


0.584 


0.81 





TABLE I. Steady-state simulation results. 
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FIG. 12. Monte Carlo simulations of coupled DP in 1 + 1 
dimensions, a): Steady-state simulations in the active phase. 
b)-d): Dynamical simulations at criticality. Time is measured 
in Monte Carlo steps. 





d= 1 


d = 2 


d = 3 


d = 4 - e 


tmax 


10 4 


10 a 


300 




5x 

&2 

5 3 


0.157(4) 

0.075(10) 

0.03(1) 


0.46(2) 
0.26(3) 
0.13(3) 


0.73(5) 
0.35(5) 
0.15(3) 


1-6/4 + 0(6") 

l/2-e/6 + 0(e 2 ) 
1/4 -O(e) 


m 
m 


0.312(6) 
3Qf9 , l 
0.47(2) 


0.20(2) 
0.56(4) 


0.10(3) 
0.43(5) 
0.75(10) 


e/12 + 0(e") 
1 /9 -I- H(V 2> I 
3/4 -O(e) 


2/zi 

2/22 
2/23 


1.26(1) 
1.25(3) 
1.23(3) 


1.10(2) 
1.12(3) 
1.10(3) 


1.03(2) 
1.04(2) 
1.03(2) 


1 + e/24 + 0(e 2 ) 



TABLE II. Dynamical simulation results. 



square spreading of the cluster vary algebraically with 
certain dynamical exponents which in turn are related 
to the exponents /3, v± 1 and m. In order to apply this 
technique to coupled DP, we prepare an initial state with 
a single A particle at the origin and perform the simu- 
lation at the multicritical point. The properties of the 
resulting cluster are analyzed separately for each particle 
species, i.e., we measure the survival probability -Pfc(t), 
the number of fc-particles (cluster mass) Nk(t), and the 
mean-square spreading from the origin Ftl(t) averaged 
over all runs that survived at level k up to time t. At the 
multicritical point, these quantities are expected to scale 
as 

P fe (t) ~ t 5 « , N k (t) ~ f fc , R 2 k (t) ~ i 2 ^ , (4.1) 

where 4 = Pk/v\\,k and z k = v\^ k /v^ k . Here, r) k is the 
so-called critical initial slip exponent [^2|,^l) which will 
be discussed below (not to be confused with the static 
correlation function Fisher exponent rf). 



The temporal variation of the quantities (4.1) mea- 
sured in a three-level coupled directed bond percolation 
model in 1 + 1 dimensions is shown in Fig. fl2"|b-d. Simi- 
lar simulations were performed in two and three spatial 
dimensions. From the slopes of the lines averaged over 
two decades we estimate the critical exponents S k , r\ k 
and Zk, which are summarized in Table Notice that 
z\, Z2, and 23 assume the same values within the nu- 
merical error. Inserting the previous estimates for (3 k , 
we can compute the scaling exponents v» k = f3 k /5 k and 
^X,fe = V\\,kl z k separately for each level fc in the hierar- 
chy (see Table [II). Within numerical error they coincide 
with the DP exponents v^, v\\, and z, as predicted by the 
field-theoretical RG calculation. 

c. General problems. The extensive simulations re- 
veal an unexpected deviation from ideal scaling at higher 
levels. As can be seen in Fig. [l^, the lines for k > 1 are in 
fact not perfectly straight but slightly bent. In order to 
illustrate these deviations, we determined the local slope 
of the survival probabilities 



5 k {t) = 



log 10 P k (2t)-\og w P k (t) 

login 2 



(4.2) 



and similarly rj k {i) in a (1 + l)-dimensional system. Af- 
ter sufficiently long times, these quantities should become 
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d= 1 


d = 2 


d = 3 


d = 4 - e 


V±,2 


1.12(4) 

1.11(15) 

0.95(25) 


0.70(4) 

0.69(15) 

0.65(15) 


0.57(4) 
0.59(8) 
0.62(9) 


l/2 + e/16 + 0(e 2 ) 


^_L,DP 


1.0968 


0.734 


0.57 




"II ,1 

"11,2 
"11.3 


1.78(6) 

1.76(25) 

1.50(40) 


1.24(6) 

1.23(17) 

1.15(30) 


1.10(8) 

1.14(15) 

1.21(15) 


l + e/12 + 0(e 2 ) 


"II, DP 


1.7338 


1.295 


1.09 





TABLE III. Derived critical exponents. 
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FIG. 13. Local slopes of the lines in Fig. [12|b-c. The slow 
drift of the slopes at higher levels indicates that the scaling 
regime is not yet reached in the present simulations. (Time 
is measured in Monte Carlo steps.) 



constant and equal to the exponents 6k and r/k, respec- 
tively. As expected, the lowest level reaches the scaling 
regime after a short time (see Fig. |l3|). At higher levels, 
however, there is a considerable drift of the local slope 
that extends over the entire temporal evolution. Simi- 
lar drifts can be observed in all other quantities which 
involve the density exponents 02,03, • ■ ■■ The deviations 
indicate that the scaling regime, especially at the third 
level, has not yet been reached. By estimating the critical 
exponents at higher levels, we therefore encounter con- 
siderable systematical errors which may even exceed the 
statistical error margins. A careful numerical analysis 
shows that the drift in the local slopes is neither related 
to finite-size effects nor to deviations from criticality. At 
present the origin of these deviations from ideal scaling 
is not yet entirely clear. It might perhaps be a signature 
of the IR-divergent diagram (c) of Fig. || (see Sect. 1IIF). 
As mentioned above, these graphs are connected with the 
appearance of additional powers of the relevant coupling 
fi, and we suspect that this drift in the scaling exponents 
signals that our first-order perturbation theory with re- 
spect to the transmutation rate becomes ultimately in- 
sufficient, and some appropriate resummation of the ex- 
pansion in /i would in fact be required. 

Furthermore, we remark that in a simulation based 
of course on a finite number of particles, ultimately 
one would expect a crossover to the decoupled situation, 
namely when the A species, whose density decays faster 
at the multicritical point, has already died out. It might 
well be possible that in a fluctuation-dominated regime 
this effect sets in much earlier, provided there emerge 
large regions which have already become depleted of the 



A particles. Thus, one explanation of the drifts visible in 
Fig. [l|| could be that this crossover region to the asymp- 
totic decoupled regime has already been reached. The 
universal exponents predicted by the field theory would 
then apply only to an intermediate scaling regime. We 
note that the field theory calculation, being based on a 
continuum description of coarse-grained particle densi- 
ties, cannot easily account for this finiteness of the par- 
ticle number. 

In the case of coupled annihilation processes (see Sec. 
VI), where similar deviations occur, the intermediate 
scaling regime can be clearly identified in numerical sim- 
ulations. In particular it is observed that the size of the 
scaling regime grows as the coupling strength decreases. 
We have also performed simulations of coupled DP with 
reduced coupling strength (probabilistic transfer of activ- 
ity to the next level) . Unfortunately, the initial crossover 
into the intermediate scaling regime grows rapidly as the 
coupling strength is reduced which makes it impossible 
to identify the boundaries of the intermediate scaling 
regime. 



B. Critical initial slip 

When a DP process starts from random initial condi- 
tions at very low density, the particles are initially sepa- 
rated by empty intervals of a certain typical size. During 
the temporal evolution these particles generate individ- 
ual clusters which are initially separated. Therefore the 
average particle density first increases as n(t) ~ f 1 — a 
phenomenon which is referred to as the critical initial slip 
of non-equilibrium systems p2fl . Later, when the growing 
clusters begin to interact with each other, the DP process 
crosses over to the usual decay n(t) ~ t~^l v ^ . Dynamical 
simulations starting from a single particle represent the 
extreme case where the critical initial slip extends over 
the entire temporal evolution. 

In ordinary DP, the critical initial slip exponent r\ is 
related to the other bulk exponents through the hyper- 
scaling relation 25 + T) = d/z |n],^lj. In the case of 
coupled DP we would therefore naively expect that the 
critical initial slip exponents r\k are related to the other 
exponents by 26k +J7fc = d/z . However, the numerical 
estimates in Tables || and III do not satisfy this scaling 
relation at higher levels k > 1. Instead they seem to 



fulfil the generalized hyperscaling relation introduced in 
Ref. [ p3| in the context of systems with many absorbing 
states, 



<5dp + 6k+Vk = d/z k 



(4.3) 



Here 6up denotes the exponent f}/v\\ of ordinary directed 
percolation. In fact, inserting the estimates of 6k and Zk 
for d = 1, Eq. (fO) predicts the values T)i = 0.314(4), 
?7 2 = 0.398(10), and n 3 = 0.443(20), which are in fair 



agreement with the estimations in Table III 



In fact, the above scaling relation (4.3) may be derived 
fairly simply, starting from an appropriate scaling form 
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FIG. 14. Susceptibility of multicritical coupled DP to an 
external field. The figure shows the stationary particle densi- 
ties rik versus the rate h for spontaneous particle creation. 



for the two point correlation function. If we initiate the 
cluster starting from a single localized seed, then the den- 
sity of B particles at a later time will have the following 
form: 



n 2 (x,t) 



\2ii 



Dt 



(4.4) 



Note that, even though this is a scaling form for the den- 
sity, it has the structure of a two point correlation func- 
tion. Roughly speaking, the prefactors in Eq. ( p~4| ) may 
be interpreted as follows: one factor of \t\ 131 comes from 
the probability that the cluster is still alive at time t, 
whilst the second factor comes from the probability that 
the point (x, t) is a mem ber of that cluster. At criticality 
we find, by integrating (4.4) over space 



~i 



(dr/X-2/30/i-ii 



f 



(Dt) 



-<t>lv\\ 



(4.5) 



wher e w e have also used one of the definitions from 
Eq. (4.1). Assuming that the scaling function / in 



the multi critical regime behaves in the same way as in 
Sec. HIE , we then end up with the scaling relation (4.2) 
for the second hierarchy level. The extension to higher 
levels works in an exactly similar manner. 



C. Susceptibility to an external field 

Coupled DP may be subjected to an external fiel d 
h > by adding a term h J d d x J dt a to the action ( |3.l| ) . 
In the particle interpretation, the field corresponds to 
a spontaneous creation of A particles at rate h during 
the temporal evolution. This means that all subsystems 
approach a fluctuating steady state, irrespective of the 
value of p. We are particularly interested in the re- 
sponse of rife to the external field at the multicritical 



point. For ordinary DP it is known that n\ ~ /i 7 , where 
7 = (3/(dv± + Uu — (3) is the susceptibility exponent. For 
coupled DP we can derive a similar relation, starting from 
the scaling form for the (steady-state) density 



n 2 



(4.6) 



where the DP exponent A can be shown to equal A 
du± + z/|| — /?!• In the limit |t| — > 0, we thus have 



n 2 



(4.7) 



Hence, using the results of Sec. HIE and generalizing to 



the fc-th level of the hierarchy, we have 



7fe = 



(3 k 



dv i 



(4.8) 



In order to verify this scaling relation, we repeat the 
steady-state simulation at the multicritical point in pres- 
ence of spontaneous creation of A particles. The results 
in 1 + 1 dimension are shown in Fig. ^4[ From the slopes 
of the lines we estimate the susceptibility exponents 

7i = 0.109(2) , 72 = 0.045(4) , 73 = 0.014(2) . (4.9) 

On the other hand, the above scaling relation yields the 
values 7i = 0.107(2), 72 = 0.051(6), and 73 = 0.018(5), 
which are in fair agreement with the simulation results. 



D. Crossover phenomena near the multicritical point 

Numerical simulations near the critical point repro- 
duce the crossover scenario predicted by the mean-field 
approximation. As an example, we consider the two 
crossovers along the dashed arrows A and A/B in the 
mean-field phase diagram of Fig. ^. To this end, we sim- 
ulate a two-level system in one spatial dimension. The 
percolation probability p± of level A is always at the crit- 
ical point, whereas level B is simulated slightly below 
and above criticality (j> 2 = p c ± 0.05). The numerical 
results are shown in Fig. [ll| As expected, for p 2 < p c 
the density of B particles first decays like n 2 (t) ~ t^ 2 /"!! 
and then crosses over to a dynamical state where the B 
particles become "slaved" to the A particles, such that 
n 2 (t) ~ ni(i) ~ t~~ ^/"ll . On the other hand, for p 2 > p c 
the B subsystem crosses over into a state with a constant 
density where the B particles become independent of the 
A particles. Thus the crossover effects are in qualita- 
tive agreement with the mean-field and RG predictions 
of Sees. ||and |HH. 



V. APPLICATIONS 

The most natural applications of coupled DP are to 
growth processes in which the layers at different heights 



22 




hi — > hi + 1 with probability q , 
or desorbed from the edge of an island (plateau) 



1000 



FIG. 15. Crossover effects near the multicritical point. 
The figure shows the particle densities rik vs time in a 
two-level system starting from a fully occupied lattice, nor- 
malized at t — 10. Level A (k = 1) is always critical, while 
level B (k — 2) is either evolving in the active, critical, or 
inactive regime. Initially the decay of B particles is the same 
in all cases. Later the system crosses over to a different be- 
havior where the B particles become independent or slaved 
to the A particles, respectively. (Time is measured in Monte 
Carlo steps.) 



represent different subsystems in the hierarchy. The dy- 
namical rules for adsorption and desorption in these mod- 
els have to be implemented in such a way that neighbor- 
ing layers are effectively coupled in one direction with- 
out feedback. The phase transition then emerges as a 
roughening transition from a smooth phase to a rough 
phase. The known examples include so-called polynu- 
clear growth models (PNG) J34| , a special class of solid- 
on-solid (SOS) models p5| ], and certain models for fungal 
growth |3(|. Another interesting realization of coupled 
DP is the spreading of activity next to the "light cone" 
in stochastic cellular automaton models with parallel up- 
date rules. 



A. Roughening transitions in SOS models 

Coupled DP was first identified in a particular SOS 
model |35f| which exhibits a roughening transition even 
in one spatial dimension. The active phase of coupled 
DP corresponds to a smooth phase where the interface is 
pinned to a spontaneously selected layer. On the other 
hand, the inactive phase of coupled DP corresponds to a 
roughening interface which propagates at finite velocity. 

The unrestricted version of the SOS model is defined 
on a one-dimensional lattice of N sites, i — 1 . . . N, with 
associated height variables hi, which may take values 
0, 1 . . . oo. The dynamical rules are defined through the 
following algorithm: At each update a site i is chosen at 
random. Then an atom is adsorbed, 



hi — ► min(hi, h i+ i) 
hi — > min(/ij, 



with probability (1 — q)/2 . 
with probability (1 — q)/2 



(5.1) 

(5.2) 
(5.3) 



Whe n the grow th rate q is low, the desorption processes 
( |5.2|) and (5.3) dominate. If all the heights are initially 
set to the same value ho, this layer will remain the bot- 
tom layer of the interface. Small islands will grow on 
top of the bottom layer, but will quickly be eliminated 
by desorption at the island edges. Thus, the interface is 
effectively anchored to its bottom layer, and the growth 
velocity, defined as the rate of increase of the minimum 
height of the interface, is zero in the thermodynamic 
limit. As q is increased, the size of islands created on 
top of the lowest layer increases. Above q c , the critical 
value of q, the islands merge and new layers are formed 
at a finite rate, giving rise to a non-zero interface velocity 
in the thermodynamic limit. 

The key feature of this model is that atoms may des- 
orb only at the edges of plateaus, i.e., at sites which have 
at least one neighbor at a lower height. In experiments 
this would correspond to a system where the binding en- 
ergy in completed layers is much larger than at the edges 
of plateaus. Furthermore, the dynamical processes at a 
given layer are independent of the processes at higher 
layers. In particular, the temporal evolution at the bot- 
tom layer is decoupled from all other processes at higher 
layers. In fact, one can show that the dynamics of the 
bottom layer can be mapped onto a one-dimensional con- 
tact process which is known to belong to the DP univer- 
sality class |35| . Identifying blank sites at the bottom 
layer as A particles, the adsorption process (5.1) may be 
interpreted as the a nnih il atio n of A particles, while the 
desorption process (5.2)-(f5~3|) corresponds to A particle 
production. Similarly, the dynamical processes at the fol- 
lowing layers may be associated with the particle species 
B,C,D,.... 

It is important to note that the state of a site in cou- 
pled DP is characterized by presence or absence of var- 
ious particle species, while sites of a growth model are 
associated with a single quantity, namely the height hi. 
To connect the two descriptions, we have to assume that 
the coupling constant a is infinite such that particles at 
level k instantaneously create particles at level k + 1. In 
this case the state of a site in coupled DP is fully char- 
acterized by the index of the lowest active level in the 
hierarchy, which then corresponds to the height of the 
interface in the growth model. Therefore, the order pa- 
rameters riA,nB, Tic-) . . . = Tii,n2,Ti3 . . . are defined by 



ho+k-l 



(5.4) 



h=ho 



that is, nk is the density of sites whose heights are less 
than h + k, where h denotes the height of the bottom 
layer. By definition, the densities obey the inequality 
nk < rik+i- 
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The above growth model is invariant under global 
shifts of the heights hi — > hi + a. This symmetry is spon- 
taneously broken in the (coupled DP) active phase where 
the system selects a particular reference height as the bot- 
tom layer. In the (coupled DP) inactive phase, however, 
the interface becomes rough and propagates at finite ve- 
locity, i.e., active DP processes subsequently enter the 
absorbing state. The growth velocity v is inversely pro- 
portional to the average survival time of the lowest-lying 
DP process, hence v ~ (q — q c ) l/ ^. Numerical simulations 
confirm that the critical behavior of the first few layers 
in the growth model is indeed the same as in coupled DP. 
In particular, the exponents 0k for the density of sites at 
the first few layers are in agreement with the numerical 
estimates in the present work. 

Alternatively, one may study the same model with the 
additional restriction \hi — hi±i\ < 1. In that case the 
layers are no longer coupled without feedback. For exam- 
ple, if a site with height hi = 1 has neighbors at heights 
hi-i = and hi+i — 2, the atom at site i cannot des- 
orb from the surface. Using the language of coupled DP, 
this means that the presence of C particles prevents the 
A particles from producing offspring. Surprisingly, the 
numerical estimates of the exponents (3k indicate that 
the critical behavior of the system is still that of coupled 
DP. Thus it seems that certain realizations of "inhibit- 
ing" feedback from higher levels to lower ones do not 
destroy the universal properties of coupled DP. Rather, 
the essential precondition for coupled DP seems to be 
the existence of a hierarchy of absorbing subspaces, i.e., 
inactive levels must not be activated by higher levels. 



B. Polynuclear growth models 

In PNG models [Q a similar scenario arises, but in 
this case the coupled DP behavior occurs at the highest 
levels of the interface. As in the previous case, the PNG 
models may be defined on a square lattice with associated 
height variables hi. The key feature of these models is 
the use of parallel updates which gives rise to a maximum 
velocity of the interface. One of the most popular PNG 
models is defined through the following dynamical rules. 
In the first half time step atoms "nucleate" stochastically 
at the surface by 



hi(t + 



1 



hi (t) + 1 , with probability p 
hi (t) , with probability 1 — p . 

(5.5) 

In the second half time step the islands grow determin- 
istically until they coalesce 



hi{t + 1) = max 
j 



(5.6) 



where j runs over the nearest neighbors of site i. Start- 
ing from a flat interface hi(0) = 0, the sites at maximal 
height hi (t) — t may be considered as active sites of a DP 



process. Obviously Eq. ( |5.5| ) turns active into inactive 
sites with probability 1 — p, w hile offspring production 
is realized by the process ( |5.6| ). Therefore, if p is large 
enough, the interface is smooth and propagates with ve- 
locity 1. Below a critical threshold, however, the density 
of active sites at the maximal height hi(t) = t vanishes, 
and the growth velocity is smaller than 1. Identifying the 
sites with hi = t as A particles, those with hi > t— 1 as B 
particles etc., the dynamical processes resemble the rules 
of coupled DP. The corresponding order parameters are 
defined by 



_^ fc-i 



(5.7) 



h=0 



Thus PNG models may be interpreted as a realization 
of coupled DP in a co-moving frame. An exact mapping 
relating PNG models and the previously discussed SOS 
models (where coupled DP resides in a fixed frame next 
to the bottom layer) was proposed in Ref. [jj5) . It should 
be emphasized that the existence of a roughening transi- 
tion in PNG models requires the use of parallel updates. 
If random-sequential updates are used, the transition is 
lost, and the interface is always rough since then there is 
no maximum velocity. 



C. Models for fungal growth 

Recently, Lopez and Jensen introduced a class of 
models for the growth of colonial organisms, such as fungi 
and bacteria. The models are motivated by recent ex- 



periments 1 37 1 with the yeast Pichia membranaefaciens 
on solidified agarose film. Depending on the concentra- 
tion of polluting metabolites, different front morphologies 
were observed. The aim of the models is to explain these 
morphological transitions on a qualitative level. 

The model for fungal growth is defined on a triangular 
(f + f )-dimensional lattice whose sites are either occu- 
pied or vacant. Growth of the colony occurs because of 
the division of individual cells, i.e., only nearest neigh- 
bors of occupied sites can become occupied. T he model 
evolves by parallel updates. To mimic realistic cells, it 
is assumed that cell division is less likely in young cells. 
To this end, the simulation keeps track of the age a,j (t) 
of occupied sites. The probability Pi(t) for a vacant site 
i to become occupied in the next time step depends on 
the total age Ai(t) — J2{ij) a jW 01 the occupied near- 
est neighbors of site i. Using the functional dependence 
Pi{t) = taah[6Ai(t)], a roughening transition was ob- 
served at 9 C = 0.183(3). Investigating clusters of sites 
growing at maximal velocity, some of the critical proper- 
ties at the transition could be related to DP |56| . It was 
argued that this roughening transition could be the es- 
sential mechanism behind the morphological transitions 
observed in experiments. 

Clearly the above fungal model and the PNG models 
are very similar in character. They both work with par- 
allel updates and exhibit a roughening transition which is 
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FIG. 16. Coupled DP in a model for fungal growth. The 
graph shows the density of the first three levels propagating 
at maximal velocity as a function of 9 — 9 C in the smooth 
phase. Power-law fits are used to estimate the exponents j3k 
(see text). The inset shows a typical configuration of the 
growing front near criticality. 



related to DP. Here we will present numerical evidence to 
show that the fungal growth model is actually a realiza- 
tion of coupled DP, in spite of complicated details such 
as the age-dependent rates for cell division and interface 
overhangs. The order parameters may be defined as 



k+l 



(5- 



where iVj(f) = 0, 1 denotes the occupation number of site 
i, and yi is the height coordinate of the i-th site. We 
have numerically measured the densities ni,ri2, and 713 
near criticality in the smooth phase (see Fig. Il6j). Our 
estimates 0i = 0.28(2), /3 2 = 0.13(2), f3 3 = 0.04(2)are in 
agreement with the numerical results of section [V . 

As in the PNG models, the existence of a roughening 
transition in the model for fungal growth requires the use 
of parallel updates. For random sequential updates there 
is no such transition and the interfaces are always rough. 
However, random sequential updates seem to be a more 
appropriate description of the experiments in Ref. |}7) , 
since realistic cells do not divide synchronously. There- 
fore it is still unclear to what extent the roughening tran- 
sition of the model in Ref. |36| is related to morphological 
transitions in realistic fungal growth. 



D. Critical behavior near the light cone in spreading 
processes with parallel dynamics 

Let us finally consider a directed bond percolation pro- 
cess on a tilted square lattice in d + 1 dimensions, which 
may be understood as a stochastic cellular automaton 
evolving by parallel updates Eaj. Starting from a single 



/ 
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FIG. 17. Realization of a unidirectionally cou- 

pled hierarchy of (1 + 1) -dimensional DP processes in a 
(2 + l)-dimensional directed bond percolation process with 
parallel updates. The figure shows the "light cone" starting 
from a single site. The subsystems A,B,C',. . . correspond to 
tilted planes as indicated by the arrows. 



active seed such a cellular automaton generates a clus- 
ter of active sites. For maximal percolation probability 
p = 1 this cluster is compact and has the shape of a pyra- 
mid. This means that all sites within the light cone (the 
surface of the pyramid) are activated. 

Apart from the usual phase transition, DP models 
with parallel updates in d > 2 spatial dimensions ex- 
hibit a second transition, where the clusters detach from 
their light cone. In the case of (2 + l)-dimcnsional di- 
rected bond percolation this transition takes place at 
p = p s ~ 0.6447 > p c . As illustrated in Fig. [l7|, the 
dynamical processes near the light cone constitute a uni- 
directionally coupled hierarchy of DP processes in d — 1 
spatial dimensions. The lowest hierarchy level A cor- 
responds to the sites in the light cone. Clearly these 
sites are decoupled from the interior of the pyramid. The 
following hierarchy levels B_,C,. . . correspond to parallel 
planes as indicated in Fig. [L7L Since activity can only per- 
colate forward in time, these planes are coupled in only 
one direction without feedback. The dynamical processes 
within each subsystem are precisely those of a directed 
bond percolation process on a tilted square lattice in d— 1 
spatial dimensions. Therefore, the numerical value of p s 
in d spatial dimensions coincides with the usual transition 
point p c in d— 1 dimensions. This explains the numerical 
value p s ~ 0.6447. 



VI. COUPLED ANNIHILATION REACTIONS 

We finally return to the question of whether new dy- 
namic universality classes can be constructed by the uni- 
directional coupling of known non-equilibrium processes. 
We have seen that in the case of linearly coupled directed 
percolation, the ensuing hierarchical structure leads to 
the emergence of multicritical behavior at a special point 
in control parameter space, described by the novel den- 
sity exponents /3j and o^. Similarly, we expect identical 
qualitative features for the closely related problem of lin- 



25 



early coupled dynamic (isotropic) percolation processes 
(albeit there one of the non-linear vertices is non-local in 
time |l4| ) . This is to be contrasted with the very general 
quadratically coupled multi-color DP processes studied 
recently by Janssen, where ordinary DP critical behavior 
is found [[15| . 

The simplest non-trivial case, however, would be to 
consider a stochastic process which is generically scale- 
invariant, i.e., where no tuning to a special critical point 
is required. An example of such a system is provided 
by the simple diffusion-limited two-particle annihilation 
reaction 



A + A 



with rate X A 



(6.1) 



The corresponding mean-field rate equation for the par- 
ticle density A reads 



dn A (x, t) 
dt 



DV 2 n A {x,t) - 2X A n A {x,t) 2 



(6.2) 



which is solved at long times t — > oo by n A (t) ~ t . 
Power counting shows that this result is expected to be 
correct for dimensions d > d c = 2. 

In low-dimensional systems, however, fluctuations and 
the emerging particle- anticorrelations become important, 
and the density decay exponent is reduced. In order to in- 
clude these fluctuations consistently, one may derive the 
following field theory from the classical master equation 



d d x / dt a{d t - DV 2 )a - X A (1 



a )a 



, (6.3) 



where we have omitted boundary terms stemming from 
the initial configuration, as well as terms related to the 
projection state (see Ref. ||). When the action (6.3) is 
expanded about the stationary solution o = 1 , the classi- 
cal field equati on f or a(x,t) yields precisely t he m ean-field 
rate equation (6.2). The entire field theory (6.3) can also 
be recast in the form of a Langevin equation for the field 
a(x,t), although this field is related to the true density 
field n A (x,t) in a rather non-trivial way [pr|,^|j3^1 . 

The structure of the field theory ( |6.3| ) is very simple, 
as no diagrams can be constructed that would renormal- 
ize the free diffusion propagator (— ico + Dq 2 )^ 1 . Fur- 
thermore, the entire perturbation series for the annihi- 
lation vertices is readily summed via a geometric series, 
or through solving the ensuing Bethe-Salpeter equation. 
Hence the scaling behavior of the density is known ex- 
actly. The final result is m 



n A (t) 




for d < 2 , 
for d — d c = 
for d > 2 . 



(6.4) 



Let us now consider a hierarchy of such annihilation 
processes, 



B + B 



with rate A# , 



(6.5) 



etc., unidirectionally coupled via the branching reaction 



A — ► A + B with rate o A b ■ 



(6.6) 



The choice of this specific coupling can be motivated as 
follows. If the A species were not to appear on the right- 
hand-side of the reaction (|6.6|), then this would constitute 
a spontaneous death process for the A particles, immedi- 
ately leading to an exponential density decay. However, 
on the lowest hierarchy level, we want to retain all the fea- 
tures of the un cou pled reactions [especially the power-law 
decay of Eq. (6^)]. Also, we want to keep the coupling 
reaction linear in the particle density n A , as in o ur earlier 
analysis of coupled DP. Thus, the reaction (6.6) feeds ad- 
ditional particles into level B, which in mean- field theory 
is described by the rate equation 



dn B (x, t) 
di 



= DX/ 2 n B {x,t) - 2X B n B {x,t) 2 
+<J AB n A (x,t) . 



(6.7) 



Obviously for long times (and consequently for low den- 
sities), the B particles are now slaved by the A species, 
and their density "adiabatically" follows n A (t), 



n B {t) 



( 



\2X l 



n A (t) 



1/2 



t 



"1/2 



(6.8) 



As is to be expected, the branching process (6.6) consid- 
erably slows down the decay on level B. Within mean- 
field theory, a straightforward generalization to higher 
hierarchy levels leads to rii(t) ~ t~ a \ with oij = 1/2 I_1 
on level i. 

In order to include fluctuation effects, we write dow n 
the action corresponding to the coupled reactions ( |6.l| ), 
(|6.5[), and (6.6), setting X A = X B = X, and a AB = a: 



S= 



d d x I dt 



a(9 t -DV>-A(l-ttV + (6.9) 



+b{d t - DV 2 )b - A(l - b z )b z + cr(l - b)aa 



Notice that the A propagator has now acquired a formal 
mass term a, as opposed to the B particles, for which we 
still the have the massless diffusion propagator (again, 
we have assumed identical diffusion constants D for both 
species). Of course, once the shifts a = 1 + a and 6=1 + 5 
are performed, this mass term disappears (as it should), 
and on the classical level the mean-field rate equation 
(|6.7|) is recovered. It is, ho weve r, convenient to work with 
the unshifted field theory ([0]), as it again has the sim- 
ple property that there are no (UV-divergent) Feynman 
diagrams that could lead to a renormalization of cither 
of the propagators. This immediately implies that the 
"mass" a is not renormalized. Thus, as opposed to the 
case of coupled DP, there is no new non-trivial scaling 
field here. If we further assume that there is no non- 
trivial contribution from the scaling function (i.e. unlike 
the case of coupled DP), th en w e may thus just insert the 
true dens ity decay results (6.4) into the mean-field rela- 
tion (6.8), in order to obtain the asymptotic decay for 
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FIG. 18. Coupled annihilation: An IR-divergent diagram 
contributing to the aaab vertex at one-loop order. 



the B particles. The result is that the decay exponents 
are halved on the second hierarchy level. The obvious 
generalization to level i is therefore 




(t-Hnt) 1 ^ 1 1 ford = d c = 2, 



for d < 2 , 

for d = d c - 
for d > 2 , 



(6.10) 



which, with the above assumption, is an exact result at 
sufficiently long times. 

However, in a simulation with finite particle numbers, 
again one would asymptotically expect a crossover to 
the decoupled scaling regime, namely when there emerge 
large regions depleted of the A species. Correspondingly, 
perhaps, one should notice that the coupled annihilation 
problem is plagued by IR-divergent diagrams which are 
very similar in nature to the coupled DP case. For exam- 
ple, to one-loop order, the newly generated aaab vertex 
for the massless shifted fields includes a diagram with 
two massless a and two massless b propagators, as de- 
picted in Fig. |l8|. This loop integral is infrared-singular 
whenever d < 6. One possible interpretation of these 
additional, apparently non-renormalizable IR singulari- 
ties, could be that they reflect an eventual non-universal 
crossover to the decoupled regime. Neither, though, can 
we exclude the possibility that ultimately a very different 
scaling regime ensues, which would have to be addressed 
by means of an effective resummation of the expansion 
with respect to the relevant coupling er. 

Numerical simulations of coupled annihilation pro- 
cesses can be performed on the same lattice as in Fig. [n|. 
The stochastic rules have to be chosen in such a way 
that a particle at site i jumps to one of the neighboring 
sites with equal probability. When two particles meet 
at the same place, they annihilate instantaneously. Sur- 
prisingly, in the limit of infinite coupling (instantaneous 
transfer of activity to the next subsystem) the resulting 
curves in a log-log plot are no t straight and do not repro- 
duce the result of Eq. flS-lOP - A more detailed analysis 
reveals that the magnitude of these deviations depends 
strongly on the coupling strength between the subsys- 
tems. To this end we replace the instantaneous transfer 



FIG. 19. The coupled annihilation process in 1 + 1 di- 
mensions: the graphs show the densities niif^jt 1 ' 2 of the first 
three levels as a function of time for different values of the 
coupling strength q. The scaling regime is marked by the two 
dashed lines (see text). Time is measured in Monte Carlo 
steps. 

of activity by a probabilistic rule, i.e., active particles 
create particles in the next subsystem at the same loca- 
tion with probability q. Clearly, q plays the role of the 
parameter a in the field theory. By v arying q we ob- 
serve that the prediction of Eq. (6.10) is only valid in 
a limited scaling regime. As q decreases, the size of the 
scaling regime grows, as illustrated in Fig. [l9|. On the 
other hand, the initial crossover into the scaling regime 
also grows with q. Similar simulations in 3 + 1 dimensions 
for maximal q suggest that these deviations still persist 
above the critical dimension although they are much less 
pronounced in that case. This supports the conjecture 
that the breakdown of the scaling regime is caused by 
IR-singular diagrams related to additional powers of er, 
which would even invalidate the simple mean-field ap- 
proach. 

Concluding this section, we note that novel critical be- 
havior does not necessarily arise in the unidirectional cou- 
pling of stochastic processes. A counterexample is given 
by the following variant of coupled annihilation, where 
we replace the reaction ( |6.6|) with 



A + A -> B + B with rate X AB ■ 



(6.11) 



The ensuing coupled diffusion-limited reaction processes 
are a special case of the more general system where the 
back reaction B + B — > A + A is present as well. The full 
system was studied by field-theoretic means in Ref. p8[ . 
Through an analysis of the coupled Bethe-Salpeter equa- 
tions for the four non-linear vertices, it was shown in |3q | 
that the A and B reactions asymptotically dec oup le, and 
each particle species decays according to Eq. (6.4). The 
physical reason for this is of course that two particles 
are required to meet in order for the coupling reaction 
Q6.ll ) to take place. Thus, this reaction competes with 
the annihilation process itself, and, in addition, as the 
daughter B particles appear on the same sites, they have 
a high probability to annihilate again immediately. This 
is somewhat related to the robustness of the DP univer- 
sality class for quadratically coupled DP processes (lSj . 
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VII. SUMMARY AND DISCUSSION 
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The simulations presented in the last few sections show 
good agreement with the predictions of the underlying 
field theory for a certain range of the parameter r for 
coupled DP. Similarly, good agreement is also found for 
a range of times t for coupled annihilation (or coupled 
DP at criticality). Nevertheless, deep into the critical 
region the simulations show a drift in the critical scal- 
ing exponents of the second and higher hierarchy lev- 
els, perhaps towards their decoupled values. It is not 
clear, however, if this drift will go all the way towards 
attaining the decoupled values of these exponents. Fur- 
thermore, the drift is more pronounced in the coupled 
annihilation model where, by decreasing the strength of 
the inter-species coupling, one can extend the range of 
the intermediate power-law behavior and delay the onset 
of the drift. 

From the field-theoretical point of view we believe 
the drift might be due to the increasing effect of the IR- 
problematic diagrams which were identified both for the 
coupled DP problem as well as for the coupled annihi- 
lation problem. These diagrams contain higher powers 
of the relevant inter-species coupling and thus are sup- 
pressed for small values of this coupling. On the other 
hand, they become more dominant for larger values of the 
transmutation rate, which, being a relevant operator, in- 
creases as one goes deeper into the critical region. This 
might perhaps render the asymptotic field theory, for 
large inter-species coupling, non-renormalizable. Note 
that the simulations for coupled annihilation show that 
the drift in the value of the exponents for the second and 
higher hierarchy levels persists even for d = 3 which is 
above the upper critical dimension d c = 2 for the first 
hierarchy level. This shows that even mean-field theory 
may not be valid at d — 3, consistent with the fact that 
there exist IR-singular diagrams diverging for any d < 6. 
Technically, a resummation of the power expansion with 
respect to \i or a would be desirable; unfortunately, a 
more satisfactory approach to this problem is not yet 
known. 

We propose the following interpretation of this sce- 
nario: eventually we expect a non-universal crossover 
into decoupled behavior. This is because in a real sys- 
tem, due to the discreteness of the number of particles, 
which is always an integer, there are likely to be large re- 
gions with no A particles at all, and there the B particles 
will behave as if they are decoupled from the A species. 
This is also true for higher hierarchy levels. Thus we 
have an interesting case in which the field theory predicts 
correctly the scaling in an intermediate universal critical 
regime, but eventually breaks down deep in the asymp- 
totic limit. There is of course the possibility that one 
might construct a meaningful field theory for the asymp- 
totic regime which will describe a crossover to a differ- 
ent critical behavior, distinct from both the intermediate 
regime and the decoupled behavior, but this seems a dif- 
ficult task and perhaps even an unlikely scenario at this 
time. 
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